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ABSTRACT: 

This  paper  describes  an  adaptation  of  the  Direct  Search  Method  which 
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straints.  This  is  sometimes  called  the  general  problem  of  mathematical 
programming.   The  algorithm  given  here  usually  performs  its  exploratory 
procedure  in  hyperplanes  approximately  tangent  to  the  constraint  hyper- 
surfaces  when  the  base  point  is  in  the  vicinity  of  such  boundaries.   There- 
fore, the  author  designates  it  the  "Tangent  Search  Method". 

Calculation  of  the  partial  derivatives  of  a  constraint  function  is 
necessary  when  it  is  violated.   However,  the  method  never  requires  evalu- 
ation of  any  derivative  of  the  objective  function.   Performance  of  Tangent 
Search  on  various  test  problems  discussed  in  this  paper  is  generally 
superior  or  similar  to  results  with  three  other  recently-published  algorithms 
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UNCLASSIFIED 


TANGENT  SEARCH  METHOD  OF  CONSTRAINED  MINIMIZATION 

Introduction 

Each  of  three  recent  articles  (l,  2,  3)  has  contributed  a  new  method 
for  solving  the  general  problem  of  minimizing  (or  optimizing)  an  arbitrary 
function  subject  to  arbitrary  constraints.   The  objective  functions  to  be 
minimized  by  these  methods  and  the  constraint  functions  considered  are  real, 
continuous,  and,  in  general,  non-linear.   However  in  practical  problems, 
for  example  "Problem  A"  of  Box  (3),  many  of  the  constraints  prove  to  be 
simple  upper  or  lower  bounds  on  independent  variables.   All  constraints 
are  expressible  as  inequalities. 

The  first  two  of  these  new  methods  can  be  classified  as  extensions  of 
the  Direct  Search  method  of  Hooke  and  Jeeves  (h) .      In  these  algorithms 
Direct  Search  (alias  "pattern  search"),  or  a  simplification  of  it,  is 
employed  until  a  constraint  is  encountered.   Then  an  alternative  or  adaptive 
procedure  is  substituted.   This  attempts  to  find  a  move  along  the  constraint 
boundary  which  will  continue  minimization  of  the  objective  function.   Trials 
are  made  by  explorations  from  the  base  point  which  always  remains  in  the 
feasible  region  where  no  constraint  is  violated.   No  trial  move  is  accepted 
as  a  new  base  unless  it  is  both  feasible  and  also  lowers  the  value  of  the 
objective  function  from  the  previous  base  value. 

In  Klingman  and  Himmelblau  (l)  the  alternative  exploratory  move  is 
called  the  "multiple-gradient  summation  technique".   Whenever  the  exploratory 
move  of  the  Direct  Search  fails  near  a  constraint  boundary  this  method  com- 
putes a  "new  successful  direction".   This  is  found  as  the  vector  sum  of  the 
normalized  gradients  of  the  contacted  constraint  and  the  objective  function. 
After  this  adaptive  move  an  attempt  is  made  to  return  to  the  Direct  Search 


algorithm  in  order  to  obviate  the  evaluation  steps  involved  in  gradient 
calculation.   The  authors  of  this  method  do  not  discuss  possible  degrading 
of  results  in  cases  when  the  gradients  can  only  be  approximated.   It  would 
seem  that  such  would  be  the  more  typical  situation  in  actual  minimization 
requirements,  such  as  Problem  A  of  Box  (3).   However,  perhaps  the  authors 
had  in  mind  recently  developed  methods  for  the  automatic  calculation  of 
derivatives,  such  as  those  reported  in  Wengert  (5)  or  Smith  (6). 

A  particularly  lucid  and  diagrammatic  explanation  of  the  reason  Direct 
Search  often  either  fails  or  becomes  inefficient  at  a  boundary  is  given  in 
Klingman  and  Himmelblau  (l).   Therefore  this  information  will  not  be  in- 
cluded in  the  present  discussion. 

In  the  Sequential  Search  method  of  Glass  and  Cooper  (2)  the  Direct 
Search  rule  is  slightly  modified.   To  keep  computationally-expensive 
exploratory  moves  to  a  minimum  their  method  always  continues  in  an  es- 
tablished successful  direction  until  failure  occurs.   That  is,  the  usual 
post-pattern  exploration  of  Direct  Search  is  omitted  except  when  the 
pattern  move  has  failed.   Consequently  it  becomes  necessary  to  build  an 
explicit  lengthening  process  into  the  pattern  move  in  order  to  accelerate 
when  an  advantageous  direction  has  been  found. 

When  the  usual  exploratory  move  fails  near  a  constraint  boundary,  the 
method  of  Glass  and  Cooper  resorts  to  an  alternative  routine  which  in 
certain  situations  must,  itself,  be  extended  to  a  further  alternative 
procedure.   They  show  that  at  such  a  base  point  the  new  required  direction 
can  be  found  as  the  solution  of  a  linear  programming  problem.   The  authors 
state  that  this  sub-problem  can  always  be  solved  by  the  Simplex  Method 
provided  a  solution  exists.   However,  the  constraint  hyper-surface  is 
often  convex,  as  defined  below.   In  this  case  the  new  direction  returns 


immediately  to  the  non-feasible  region.   Therefore  an  "extension"  to  the 
alternative  move  provides  systematic  rotation  of  the  indicated  direction 
until  it  clears  the  constraint.   After  an  alternative  move  an  attempt  is 
made  to  return  to  the  modified  Direct  Search  algorithm,  since  there  is 
obviously  considerable  expenditure  of  computational  effort  in  carrying 
out  the  procedures  outlined  above. 

Unfortunately  only  two  very  simple  examples,  each  with  only  two  in- 
dependent variables,  were  included  in  Cooper  and  Glass  (2).   Therefore 
the  present  author  is  not  convinced  that  it  offers  an  efficient  method  of 
solving  practical  problems. 

The  Complex  Method  of  M.  J.  Box  (3)  falls  into  a  different  category. 
It  is  a  novel  and  elegant  adaptation  of  the  Simplex  Method  of  linear 
programming,  and,  therefore,  is  not  conceptually  allied  to  the  pattern 
search  basis  of  the  two  methods  outlined  above.   Starting  with  any  one 
feasible  first  point  in  n-dimensional  space  a  "complex"  of  2n  vertices  is 
constructed  by  selecting  random  points  from  the  feasible  region.   Then  one 
simple  computational  loop  is  employed.   These  instructions  find  the  current 
worst  vertex,  that  is,  the  vertex  with  largest  corresponding  value  for  the 
objective  function,  and  replace  it  by  its  oVer-ref lection  through  the 
centroid  of  all  other  vertices.   If  the  vertex  to  be  replaced  is  considered 
as  a  vector  in  n-space,  its  overf lection  is  opposite  in  direction,  increased 
in  length,  and  collinear  with  the  old  vertex  and  the  centroid  of  the  other 
vertices.   When  the  overf lection  is  not  feasible  or  remains  worst,  it  is 
displaced  half  way  toward  the  centroid.   Constraints  are  classified  as 
explicit  or  implicit.   The  former  are  simple  upper  or  lower  bounds  on  in- 
dependent variables.   The  latter  are  treated  as  upper  or  lower  bounds  on 
arbitrary  functions  of  the  independent  variables.   This  distinction  is 


used  profitably  in  the  Box  algorithm. 

Eventually  all  vertices  converge  to  a  point  which  is  taken  as  the  so- 
lution.  Preliminary  experience  seems  to  show  that  the  Complex  Method  is 
very  dependable  but  not  completely  infallible.   (See  problems  3,  5,  7  and 
9  in  the  numerical  examples  below. ) 

Basic  Concepts  of  the  Tangent  Search  Method. 

The  Tangent  Search  Method  which  is  described  in  the  present  paper  falls 
into  the  pattern  search  family  of  suggested  solutions  to  the  problem  of 
constrained  minimization.   The  modified  pattern  search  of  Glass  and  Cooper 
is  used  except  that  the  test  for  constraint  violations  is  postponed  until 
the  exploratory  move  is  completed  for  all  independent  variables.   Whenever 
a  trial  move  passes  a  boundary  the  partial  derivatives  of  the  violated  con- 
straint function  are  approximated  at  the  current  feasible  base  point.   This 
is  performed  by  the  most  elementary  difference  method  possible.   Every 
effort  is  made  to  eliminate  unnecessary  use  of  the  auxiliary  procedure  which 
evaluates  constraint  functions. 

Then  an  exploratory  move  is  made  in  the  hyperplane  which  is  approximately 
tangent  to  the  constraint  hyper-surface.  (The  hyperplane  is,  in  general, 
slightly  removed  into  the  feasible  region  so  that  it  contains  the  current 
base. )   This  substitute  exploratory  move  seemed  to  the  present  author  a 
more  obvious  adaptation  of  the  Hooke  and  Jeeves  algorithm  than  either  of 
the  two  described  above.   If  it  is  successful  it  is  called  a  "tangent  move". 
Economy  over  the  method  of  Klingman  and  Himmelblau  is  achieved  since  evalu- 
ation of  the  gradient  of  the  objective  function  is  completely  avoided.   The 
difference  technique  used  to  approximate  the  partial  derivatives  of  the 
constraint  function  will  be  described  in  a  separate  section  before  the 
numerical  examples  are  considered. 


An  additional  feature  of  the  algorithm  here  described  is  the  allowance 
made  for  small  perturbations  of  the  base  point  when  the  tangent  exploration 
procedure  fails  several  times  consecutively.   Under  some  circumstances  it 
seems  better  to  transfer  provisionally  the  base  point  to  a  higher  feasible 
position  on  the  objective  function  than  to  continue  reducing  step  sizes. 
This  type  of  move,  called  "jump  move"  below,  helps  to  extract  the  base 
point  from  difficult  slots  and  corners  in  the  feasible  region  caused  by 
the  multiplicity  of  constraints. 

Finally,  it  is  obvious  that  the  directions  for  exploration  in  the 
tangent  hyperplace  can  be  chosen  in  any  number  of  ways.   The  attempt  of 
this  method  is  to  preserve  the  basic  logic  of  the  usual  exploratory  move, 
but  to  allow  a  variety  of  possibilities  within  that  outline.   Previous 
success  in  finding  a  good  "tangent  move"  for  a  constraint  dictates  the 
mode  of  the  tangential  exploration  first  tried  if  that  constraint  is 
violated  again. 

Comparisons  of  the  performance  of  the  proposed  method  and  those  pre- 
viously mentioned  will  be  given  in  the  penultimate  sections  of  this  paper. 
The  immediately  following  sections  give  an  informal  mathematical  definition 
of  the  problem  of  constrained  minimization  and  then  a  detailed  exposition 
of  the  Tangent  Search  Method. 

The  Problem  of  Constrained  Minimization. 

Let  X  =  ( x, ,  Xp,  ...,  x_)  be  a  vector  in  n-  dimensional  space.   Let 
f  (X)  and  g.  (X),  (j  =  1,  2,  ...,  m) ,  be  continuous  functions  which  can 

J 

be  evaluated  for  any  X  within  the  vector  space. 


Any  X  for  which 

g   (X)  >  0,   (j  =  1,  2,  ...,  m), 

J 

will  be  called  "feasible". 

Also,  any  constraint  hyper-surface, 
gk  (X)  =  0, 
will  be  called  "convex"  in  a  region  of  interest  if  every  point  of  every 
hyperplane  tangent  to  that  constraint  hyper-surface  in  the  area  of  interest 
is  not  feasible. 

The  general  problem  of  constrained  minimization  is  to  determine  the 

vector,  X  ,  such  that   f(Xs)  is  less  than  all  other  f(X),  subject  to  the 

constraints,  g.  (X  )  ^0,  for  all  j  =  1,  2,  ...,  m. 
J    s 

If  for  any  point,  X, 
g  (X)  <  0,  then 

J 

the  jth  constraint  is  "violated"  at  X. 

For  the  problems  of  constrained  minimization  considered  here, 
gk(X  )  =  0,  where  1  £  k  4  m. 
That  is,  typically,  the  solution  will  lie  on  one  or  more  constraint  boundaries. 

An  important  requirement  of  the  Tangent  Search  Method  is  that  the  first 

partial  derivatives  of  g.(X)  with  respect  to  each  X.,  (i  =  1,  2,  ...,  n) 

J  ■*- 

should  be  continuous  and  computable,  at  least  in  the  vicinity  of  any 
violation.   Special  provision  is  made  for  cases  where  such  partial  deriva- 
tives are  zero. 

Details  of  the  Presentation:  Chart  Symbols  and  Notes 

The  discussion  below  will  assume  familiarity  with  the  method  of 
Direct  Search  as  given  by  Hooke  and  Jeeves  (U).   It  is  also  outlined  in 
Glass  and-  Cooper  (2)  and  Klingman  and  Himmelblau  (l).   The  Tangent  Search 


algorithm  will  be  communicated  by  means  of  detailed  flow  diagrams  similar 
to  those  of  Hooke  and  Jeeves.   Descriptive  flow  diagrams  will  be  omitted. 
However  each  box  of  the  detailed  flow  diagrams  will  be  numbered  within 
parentheses.   Corresponding  descriptive  titles  (underlined)  and  additional 
explanatory  notes  will  be  given  in  the  text  following  each  chart. 

Notation  will  differ  from  that  of  previous  papers  in  order  to  avoid 
Greek  symbols  which  are  inconvenient  on  a  typewriter.   Flow  chart  con- 
ventions will  be  more  standard: 

1.  The  diamond  will  represent  a  decision  function. 

2.  The  circle  will  represent  a  connector.   Off-page  connector  numbers 
will  be  preceded  by  the  referenced  chart  number(s);  for  example, 
connector  5  of  chart  3  would  be  denoted  "3/5"- 

3.  The  square  will  represent  starting/stopping  points  in  the  main 
procedure  and  entry/exit  points  for  auxiliary  procedures. 

k.      Processing  functions  are  represented  by  rectangles. 

5.  Input/output  functions  are  shown  as  rectangles  with  rounded 
corners.   Only  the  most  basic  of  these  are  given. 

6.  When  an  operation  involves  a  matrix,  indicated  by  capitalization, 
it  should  be  understood  that  it  is  carried  out  on  all  elements. 
Otherwise  subscripts  indicate  the  elements  referenced. 

7.  The  "equals"  sign  is  used  as  the  replacement  symbol;  the  evaluated 
right  hand  side  replaces  the  former  value  (if  any)  of  the  variable 
denoted  on  the  left  hand  side. 

Auxiliary  Procedures 

Two  external  auxiliary  procedures  are  used  in  various  parts  of  the  flow 
diagrams.   The  actual  forms  of  these  vary  according  to  the  objective  function 


to  be  minimized  and  the  constraints  imposed  on  the  solution.   The  notation 
and  purpose  of  each  are  given  below: 

1.  F(x)  is  a  function  procedure  which,  when  called,  produces  a  value 
of  the  objective  function  at  any  point,  X,  whether  or  not  X  is 
feasible. 

2.  C(X)  is  a  subroutine  which,  when  called,  produces  a  column  matrix 

of  m  elements.   These  are  the  values  at  the  point  X  of  the  constraint 

functions,  g  ,  for  j  =  1,  2,  ...,  m. 
J 

The  eight  internal  auxiliary  procedures  are  not  changed  from  problem  to 
problem.   They  are  separated  from  the  main  procedure  in  the  discussion  below 
either  because  they  are  called  from  more  than  one  site  or  because  logical 
delineation  of  the  function  performed  may  clarify  the  presentation.   The 
notational  representation  for  the  internal  auxiliaries  is  adapted  from 
Hooke  and  Jeeves  (k) .   A  colon  separates  the  auxiliary  name-symbol  from 
its  parameter  list.   The  list  is  divided  between  input  and  output  parameters 
by  a  semicolon.   As  an  example,  the  notation 

"Q  :  x,  X;  y,  Y", 
would  indicate  that  the  internal  auxiliary,  Q,  using  current  values  of  the 
variables,  x  (a  scalar)  and  X  (a  vector  or  matrix),  causes  replacement  of 
the  values  of  the  variables,  y  (a  scalar)  and  Y  (a  vector).   If  an  argument 
is  both  an  input  and  output  parameter,  it  will  appear  only  in  the  latter 
list. 

The  details  of  the  internal  auxiliary  procedures  E,  E  ,  F  ,  H,  T, 
and  V  are  given  in  flow  diagrams  of  later  sections. 

The  internal  auxiliary  procedure,  M,  is  simply  the  usual  minimum-value 
function:  e_.g_.  , 
"M:X;  x" 


indicates  that  the  algebraically  smallest  element  of  X  replaces  the  present 
value  (if  any)  of  x. 

The  internal  auxiliary  procedure,  S,  is  the  usual  sign  function;  e_.g_.  , 
"S  :  x,  y;  z" 
indicates  that  the  value  of  z  is  replaced  by  the  absolute  value  of  x, 
prefixed  with  the  sign  of  y. 

Matrix  Variable  Symbols 

A   Values  of  the  m  constraint  functions  at  current  best  feasible 

move  produced  by  auxiliary  procedure  T. 
B   Values  of  the  m  constraint  functions  at  base  point. 
D   Current  exploration  step  sizes  for  the  n  variables. 
G   Values  of  the  m  constraint  functions  at  trial  move. 
I   Matrix  of  indices  containing  one  column  of  n  variable  numbers 

for  each  of  the  v  constraints  presently  violated. 
J   Auxiliary  column  matrix,  defined  in  the  same  manner  as  K. 
K   Indices  of  the  v  currently-violated  constraints. 
L   Lower  limits  allowed  on  the  corresponding  n  values  of  D.   When 

Dj_<  L.  ,  for  all  i  =  1,  2,  ...,  n,  the  search  is  terminated. 
N   Table  containing  the  current  coupling  mode  for  each  of  the  m 

constraints . 

P   Matrix  of  v  columns,  each  containing  n  first  partial  derivatives 

i 
of  a  violated  constraint  function  with  respect  to  the  independent 

variables.   These  are  arranged  in  ascending'order  of  absolute  value, 
U   Upper  limits  allowed  on  the  corresponding  n  values  of  D. 
W   Base  point  prior  to  last  jump  move,  containing  n  values  of  the 

variables . 


10 


X   Trial  move  (n  elements).   Also,  the  vector  of  independent 
variables  of  n-space  in  mathematical  discussion. 

Y   Current  base  point  (n  elements). 

Z   In  the  main  procedure,  the  previous  base  point.   It  is  also  used 
as  a  temporary  array  by  auxiliary  procedures,  H  and  T. 

Scalar  Variable  Symbols 

a   Index  of  primary  variable  of  any  component  of  a  tangent  explora- 
tory move. 

b   Index  of  secondary  variable  in  any  component  of  a  tangent  explora- 
tory move. 

c   Failure  counter;  controls  application  of  jump  move  option. 

d   Constant  used  in  step  size  control. 

d^  Decrement  of  the  secondary  variable  in  any  component  of  a  tangent 
exploratory  move. 

e   Temporary  indicator,  used  variously. 

f   Value  of  objective  function  after  trial  move.   Also  denotes  the 
objective  function  in  the  mathematical  discussion. 

f   Value  of  objective  function,  corresponding  to  tentative  jump  move. 

fw  Value  of  objective  function  at  base  prior  to  last  jump  move. 

fx  Function  value  used  temporarily  by  conventional  exploratory  move. 

f   Value  of  objective  function  at  current  base  point. 

g   Coupling  constant  for  active  variables  in  auxiliary  procedure  T. 

g.   The  j-th  constraint  function  (in  mathematical  discussion). 
J 

h   Coupling  constant  for  passive  variables  in  auxiliary  procedure  T. 
i,  j,  k  Utility  indices  used  variously, 
m   Number  of  constraints, 
n   Number  of  independent  variables. 

11 


p.   First  step  extension  factor. 

p   Second  step  extension  factor. 
2 

p   Partial  derivative  of  constraint  function  with  respect  to  primary 
a 

variable  of  any  component  of  auxiliary  procedure  T. 

p   Partial  derivative  with  respect  to  secondary  variable.   (See 
b 

p  above. ) 
a 

q  Temporary  storage  and  utility  index. 

r  Step  reduction  factor. 

s  Utility  switch,  used  variously. 

t  Feasible  move  indicator. 

u  Auxiliary  counter,  used  in  the  same  manner  as  v. 

v  Number  of  constraints  currently  violated. 

w  Coupling  variable  index,  used  in  auxiliary  procedure  T. 

x  Variable  increment  used  in  auxiliary  procedure  H. 

x  ,  x   Primary  and  secondary,  respectively,  trial  increments  at  any 
a   b 

component  of  tangent  exploratory  move. 

x   The  i-th  variable  in  n-space  in  mathematical  discussion, 
i 

y  ,  y   Previous  values  of  primary  and  secondary  variables  before  trial 
a  tq 

component . 
z   Variable  increment  used  in  auxiliary  procedure  H. 

Details  of  the  Method;  Initialization  and  Search  for  Feasible  Start 

The  detailed  flow  diagram  for  these  phases  of  a  solution  is  given  in 

Chart  1.   At  times  it  is  a  considerable  convenience  to  require  that  the 

algorithm,  itself,  find  a  feasible  initial  point  before  the  actual  search 

for  constrained  minimization  is  initiated.   This  is  easily  accomplished 

by  using  a  variation  of  the  usual  unconstrained  Direct  Search  to  minimize 

the  function  F   :  X,  G;  f  .   This  internal  auxiliary  procedure,  which  is 
o  x 
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Chart   1:      Initialization  and  Search  for  a  Feasible  Start 
Main  Procedure,   Part   I 
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shown  in  Chart  3,  simply  calculates  the  sum  of  the  squares  of  all  constraint 
functions  which  are  negative  at  point  X.   As  soon  as  f  "becomes  zero,  the 
actual  solution  can  begin. 

Descriptive  titles  for  each  numbered  block  of  Chart  1  and  additional 
explanatory  notes  are  given  below. 

(1)  Input  starting  point ,  computational  constants ,  initial  and  minimal 
step  sizes. 

Optimal  choice  of  r,  the  step  reduction  factor,  has  not  been 
thoroughly  investigated.   The  example  problems  discussed  in  a  later  section 
were  run  with  r  =  .25,  .  375»  and  .5-   These  three  values  seemed  to  produce 
approximately-equivalent  satisfactory  results  for  various  problems  in 
unconstrained  minimization  previously  solved  by  the  author  in  trials  of 
the  conventional  Direct  Search  algorithm  on  the  Control  Data  l6oU  computer. 
However,  as  will  be  shown  in  the  discussion  of  numerical  examples  below- 
significant  differences  in  performance  developed  in  some  problems  of 
constrained  minimization.   According  to  Hooke  and  Jeeves  (U)  the  simple 
Direct  Search  algorithm  in  unconstrained  problems  is  not  sensitive  to  this 
choice,  as  long  as  r  <  1. 

(2)  Initialize  computational  constants. 

The  pattern  and  exploration  step  size  extension  factors  are  both 
1  +  r;  however,  the  exploration  steps  are  limited  to  a  maximum  size  of 
four  times  the  original  values. 

(3)  Evaluate  objective  function  and  constraint  functions  at  starting 
point;  check  for  violations. 

The  coupling  modes  are  initialized  at  zero  for  all  constraints.   The 
use  of  these  quantities  is  explained  in  a  subsequent  section.   Details  of 
the  internal  auxiliary  procedure  V,  are  given  in  Chart  2.   Briefly  stated, 
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this  procedure  examines  the  constraint  functions,  counts  the  number  of 
violations,  and  stores  the  index  number  of  each  violated  constraint  in  K. 
( h )   Is_  any  constraint  violated? 

If  not,  exit  to  the  actual  constrained  search  is  made.   If  so,  a 
feasible  starting  point  must  first  be  found. 

(5)  Evaluate  special  function  at  origin. 

Details  of  this  internal  function,  F  ,  are  given  in  Chart  3. 

(6)  Starting  at  base  point,  perform  exploratory  move,  E  ,  using  the 

special  function,  F  . 

— *■ o 

A  detailed  flow  diagram  of  the  Hooke's  and  Jeeves'  exploratory  pro- 
cedure, E,  is  given  in  Chart  h.      The  auxiliary,  EQ,  is  identical  to  E, 
except  that  all  appearances  of  "f  =  F(X)"  are  to  be  replaced  by 
"fx  =  F0(X)". 

( 7 )  Is.  trial  value  of  special  function  zero? 

If  so,  a  feasible  starting  point  has  been  found. 

( 8 )  Is_  trial  value  of  special  function  below  the  base  value ? 
If  so,  a  pattern  move  can  be  attempted. 

( 9 )  Are  all  step  sizes  below  corresponding  minimal  criteria? 

(10)  Stop 

If  control  arrives  at  Item  (10),  Direct  Search  has  failed  to  locate 
any  feasible  starting  point  for  the  actual  search.   Perhaps  the  constraints, 
as  defined  in  external  auxiliary  procedure  C(X),  are  inconsistent,  or  the 
elements  of  L  are  too  large. 

(11)  Reduce  exploratory  step  sizes. 
Another  exploratory  move  follows . 
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Chart  2:   Internal  Auxiliary  Procedure  V:  G;v,K. 

This  procedure  places  a  count  of  currently-violated 
constraints  at  v  and  tabulates  their  indices  in  col- 
umn vector  K. 
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Chart   3:      Internal  Auxiliary  Procedure  F   :    X,G;f    . 
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Chart  4j   Internal  Auxiliary  Procedure  E:  Y,D;  f ,X. 

This  is  the  conventional  exploratory  algorithm  of 
Hooke  and  Jeeves. 
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(12)  Set  new  base  point;  make  extended  pattern  move;  evaluate 
special  function. 

(13)  Is.  trial  value  of  special  function  zero? 

If  so,  pattern  move  arrived  at  feasible  start. 
( Ik )   Is_  trial  value  of  special  function  below  base  value? 
If  not,  an  exploratory  move  follows. 

( 15 )  Increase  exploratory  step  sizes,  subject  to  an  upper  limit . 
A  successful  pattern  move  causes  these  steps  to  be  extended. 

Another  extended  pattern  move  follows. 

( 16 )  Evaluate  constraint  functions  at  feasible  start ;  set  this  point 

as  first  base;  restore  original  step  sizes;  evaluate  objective 

function. 

Initialization  for  the  constrained  search  is  accomplished  by  these 

\ 
steps . 

The  preceding  three  charts  show  the  details  of  three  important 
internal  auxiliary  procedures.   They  do  not  seem  sufficiently  complicated 
to  warrant  inclusion  of  descriptive  titles  and  additional  notes. 
Modified  Direct  Search 

Chart  5  can  be  conveniently  considered  as  representing  the  steps  of 
the  algorithm  which  are  followed  when  the  base  point  is  not  in  the 
vicinity  of  any  constraint  boundary.   The  method  shown  is  substantially 
Direct  Search  (U)  with  some  modifications  very  similar  to  those  proposed 
by  Glass  and  Cooper  (2). 

It  should  be  recalled  that  Direct  Search  assumes  an  initial  feasible 
X,  and  the  corresponding  f  is  first  computed.   In  order  to  explore, 
each  X.  is  increased  by  its  D.  and  if  the  corresponding  f  is  below 
f,  the  trial  X.  is  adopted.   If  not,  the  effect  of  decreasing  X^   by 
D.  is  tried.   If  both  +_  D.  components  fail  to  lower  the  objective 
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function,  X.  is  restored  to  the  base  value.   The  exploratory  move  tries 
1 

all  variables  in  this  manner,  and  if  the  complete  move  produces  an  im- 
proved f,  the  vector  between  the  base  point  and  the  successful  exploratory 
move  provides  the  direction  for  a  trial  pattern  move  and  its  reference 
length.   The  first  pattern  step  in  that  direction  is  (l  +  r)  times  this 
reference  length.   Accelerating  moves  are  made  in  the  same  direction  until 
some  new  trial  value  of  f  is  greater  than  the  base  value,  f  . 

When  a  pattern  or  exploratory  move  fails,  every  element  of  D  is  reduced 
by  a  multiplicity  factor,  r  4.  1.   Eventually  either  a  new  successful 
pattern  is  found  or  else  the  lower  limits,  L,  on  all  D  are  passed.   When 
this  happens  the  calculation  is  terminated. 

The  method  of  pattern  move  acceleration  differs  from  Hooke  and  Jeeves. 
The  present  author  believes  that  this  variation  usually  reduces  the 
required  number  of  exploratory  moves  which  are  costly  in  terms  of  objective 
function  evaluations.   This  conclusion  was  also  reached  by  Glass  and 
Cooper  (2).   Furthermore,  it  seems  to  follow  logically  from  the  following 
observation  on  page  2l6  of  Hooke  and  Jeeves  (h) .      "Typically  a  pattern 
once  established  will,  through  continuous  modification,  grow  until  the 
length  of  the  pattern  move  is  10  to  100  times  the  basic  step  size." 
Any  successful  pattern  move  leads  to  a  longer  one  in  exactly  the  same 
direction;  it  also  causes  a  corresponding  increase  in  every  element  of  D. 
However,  upper  limits  on  the  value  of  D  are  imposed.   These  processes  are 
also  not  part  of  the  original  Direct  Search.   They  seem  justified  on  the 
basis  of  the  favorable  results  achieved. 

The  fundamental  variations  from  the  method  of  Hooke  and  Jeeves  are, 
of  course,  inclusion  of  tests  for  constraint  violations  and  the  alternate 
exploratory  procedure  employed  when  such  occur.   Tests  for  violations  are 
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made  after  each  otherwise-successful  exploratory  and  pattern  move.   When- 
ever one  or  more  boundaries  have  been  crossed,  exit  is  made  to  the  tangent 
exploratory  procedure.   This  will  be  described  in  the  later  sections. 

Finally,  the  jump  move  is  proposed  as  a  wholly  new  tactic  in  finding 
difficult  solutions.   If,  at  one  or  more  constraint  boundaries,  the  tangent 
exploratory  procedure  repeatedly  fails  to  establish  a  successful  direction 
(tangent  move),  a  jump  move  is  made  to  the  best  feasible  trial  of  all 
components  developed.   That  is,  the  base  point  is  shifted  to  that  com- 
ponent of  the  v  tangent  exploratory  moves5  just  attempted  by  auxiliary 
procedure  T,  which  is  feasible  and  has  the  lowest  corresponding  value  of 
the  objective  function. 

The  solution  then  proceeds  with  an  initial  reapplication  of  the  tangent 
exploratory  procedure.   If  repeated  failures  at  one  or  more  constraints 
recur,  another  jump  move  is  made,  provided  that  the  intervening  solution 
steps  have  produced  a  base  point  with  function  value  lower  than  the  value 
prior  to  the  first  jump  move.   If  this  has  not  happened,  the  subsequent 
potential  jump  move  is  bypassed.   Instead,  the  base  is  returned  to  its 
position  before  the  first  jump  move  and  a  conventional  exploratory  move 
is  attempted.   More  precise  information  on  this  feature  will  be  given  in 
the  following  charts  and  notes . 
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Chart  5:  Modified  Direct  Search  with  Jump  Move. 
Main  Procedure,  Part  II 
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Notes  for  Chart  5 

( 1 )  Set  failure  counter  to  zero;  starting  from  base  point,  perform 
exploratory  move. 

The  details  of  this  conventional  exploration  were  shown  in  Chart  U, 
above . 

(2)  I_s_  trial  objective  function  below  value  at  base? 

If  so,  constraints  must  next  be  checked  for  any  possible  violation. 

( 3)  Set  failure  counter  to  n  +  1. 
(h)      Reduce  step  sizes. 

This  is  also  the  point  of  re-entry  into  the  main  procedure  after 
failure  of  the  tangent  explore  move. 

( 5 )  Are  all  step  sizes  less  than  lowest  tolerated  limit? 
If  so,  the  calculation  is  terminated. 

(6)  Output  values  of  variables  and  objective  function  at  minimum. 

The  accuracy  obtained,  provided  the  method  is  successful,  depends 
upon  the  specified  values  of  the  L  column  matrix. 

( 7 )  Set  and  test  switch  s_.   If_  s  £  0,  try  tangent  exploration;  if 
s  =  0,  test  for  existence  of  &   feasible  move;  if  s>0,  try 
conventional  exploration. 

Only  after  n  consecutive  applications  of  the  tangent  exploration 
algorithm  fail  to  produce  a  satisfactory  move  does  a  jump  move  become 
possible.   Switch  s  will  be  greater  than  zero  only  immediately  after 
failure  of  a  conventional  exploration. 

( 8 )  Has  a   feasible  jump  move  been  stored  by  T? 

If  the  previous  application  of  auxiliary  procedure  T  produced  no 
feasible  trial,  it  is  re-entered. 
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(9)  ls_   the  objective  function  at  the  present  base  below  its  value 
before  the  most  recent  jump  move? 
This  is  a  necessary  condition  for  any  subsequent  jump  move. 

( 10 )  If_  not,  restore  previous  base  and  try  conventional  exploratory  move. 

(11)  If_  so,  perform  jump  move;  then  go  to  recompute  partial  derivatives 
before  tangent  exploration. 

The  present  base  values  are  saved,  so  that  the  test  at  Item  (9)  can 
be  made  the  next  time  a  jump  move  may  be  indicated.   The  best  feasible 
trial  from  the  most  recent  application  of  procedure  T  becomes  the  new  base, 
and  the  present  base  is  treated  as  a  failed  trial. 

( 12 )  Compute  constraint  functions  and  check  for  violations. 

( 13 )  _Is_  there  any  violation? 

If  so,  control  is  transferred  to  the  preamble  to  the  tangent 
exploratory  procedure. 
( lU )   Perform  extended  pattern  move. 

The  failure  counter  is  set  to  unity  as  an  indicator;  present  values 
of  the  constraint  functions  are  saved  for  possible  use  by  auxiliary  pro- 
cedure H. 

( 15 )  Is_  the  trial  value  of  the  objective  function  below  base  value? 
If  not,  pattern  move  has  failed  and  exploration  is  indicated. 

( 16 )  Extend  exploration  step  sizes  . 

Whenever  a  pattern  move  succeeds  in  lowering  the  objective  function 
the  values  of  D  are  increased,  subject  to  a  limit. 

The  Exploratory  Move  in  the  Tangent  Hyperplane;  Basic  Concepts. 

Whenever  at  some  trial  point,  X,  a  conventional  exploratory  or  pattern 
move  finds  a  value  of  the  objective  function  lower  than  the  base  value, 
that  trial  move  is  immediately  checked  to  determine  whether  or  not  it  is 
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feasible.   If  not,  control  is  transferred  to  the  preamble  of  the  tangent 
exploratory  procedure,  shown  in  Chart  6.   This  section  will  provide  a 
general  outline  of  tangent  exploration,  including  ancillary  logic  such 
as  the  preamble  and  the  calculation  of  partial  derivatives  of  the  con- 
straint hyper-surfaces.   Later  sections  will  explain  each  sub-procedure 
in  greater  detail. 

The  internal  auxiliary  procedure  H:v,  B,  G,  K,  X,  Y;  I,  P  approximates 
the  partial  derivative  of  each  constraint  function  currently  violated  at 
the  last  feasible  base,  i.e_.  ,  the  quantities; 
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(j  =  1,  2,  ...,  v),  (i  =  1,  2,  ...,  n) 


It  also  generates  the  matrix  I  which,  after  exit  from  H,  contains  a 
column  of  indices  for  each  violated  constraint.   Each  of  these  v  columns 
contains  n  variable  numbers  arranged  in  order  of  ascending  absolute  value 
of  the  corresponding  first  partial  derivatives.   If  the  tangent  exploratory 
move  fails,  auxiliary  H  is  not  re-entered  on  subsequent  attempts  to  find 
a  tangent  move  from  that  base.   If,  however,  after  n  such  failures  a  jump 
is  made,  H  is  used  to  re-estimate  the  appropriate  derivatives  before 
procedure  T  is  re-entered. 

If  Y  is  the  current  base  point  in  feasible  n-space  and  the  j-th  con- 
straint hyper-surface  is  crossed  due  to  trial  move  X,  the  first  essential 
concept  of  the  Tangent  Search  Method  is  that  the  next  exploration  be  made, 
if  possible,  in  the  hyperplane, 
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where  x.  is  the  i-th  independent  variable  of  n-space. 

Note  that  this  hyperplane  is  only  approximately  tangent  to  the  hyper- 
surface,  g.(x)  =  0,  since  it  includes  the  point  Y  which,  in  general,  is  dis- 
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placed  some  amount  into  feasible  space.   However,  the  discussion  following 
will  usually  refer  to  it  as  the  "tangent  hyperplane".   It  would  be  so 
strictly  only  if  the  base  lay  directly  upon  the  constraint  boundary. 

The  second  basic  idea  of  the  tangent  exploratory  scheme  is  that  each 
variable,  in  turn,  be  incremented,  while  holding  all  others  but  one  con- 
stant.  That  one  is  adjusted  by  the  amount  necessary  to  confine  this  trial 
component  of  the  move  to  the  tangent  hyperplane.   The  a-th  variable,  in- 
cremented by  the  corresponding  D  ,  will  be  called  the  "primary"  variable 

a, 

of  a  trial  component  of  a  tangent  exploratory  move.   The  b-th  variable, 

which  is  that  one  adjusted  to  keep  the  trial  in  the  hyperplane,  will  be 

called  the  "secondary"  variable.   Each  such  trial  component,  therefore, 

requires  the  two  replacement  operations: 

X  =  X  +  D 
a    a    a 


Xv  =  X 

b     b 


^gj 

^xa 

Y 

^gj 

The  primary  index  simply  cycles  through  the  first  (n  -  l)  elements 
of  the  j-th  column  of  the  current  I  matrix,  i.e.,  a  =  I,  . ,  (k=l,  2,  ..., 
n-l).   The  corresponding  b  indices  are  selected  by  a  rule  which  will  be 
discussed  in  the  section  below  treating  the  preamble  to  tangent  exploration. 

As  in  the  conventional  exploratory  move,  the  objective  function  is 

evaluated  for  any  trial  component,  the  k-th,  for  example.   If  it  is  lower 

than  any  previous  trial  of  the  present  move  and  if  no  constraint  is  violated, 

that  component  is  adopted  and  the  variable  indexed  by  I,      .   becomes 

(k  +  l)j 
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primary. 

Otherwise  the  reflection  of  that  trial  component  is  tested,  i_.e_.  , 
the  component  generated  by  the  replacements, 

and 

D  . 
Y  a 

On  the  first  component  this  reflection  is  collinear  with  the  original  trial 
and  the  base  point.   On  subsequent  components  it  is  collinear  with  the 
previous  trial  and  the  best  move  yet  developed  during  the  current  ex- 
ploration.  If  this  trial  component  is  feasible  and  produces  a  function 
value  lower  than  the  previous  trial,  it  is  adopted.   Otherwise  both  the 
a-th  and  b-th  variables  are  restored  to  their  values  before  the  former 
became  primary,  and  the  next  primary  variable  is  substituted.   This 
process  continues  until  (n  -  l)  primary  variables  have  been  selected. 

Preamble  to  Tangent  Exploratory  Move. 

Introductory  steps  before  entrance  into  and  testing  procedures  after 
exit  from  procedure  T  are  shown  in  Chart  6. 

If  the  violation  which  results  in  tangent  exploration  was  caused  by 
conventional  exploration,  step  sizes  are  reduced  in  the  preamble.   How- 
ever, if  an  attempted  pattern  move  caused  the  violation,  step  sizes  are 
unchanged.   This  policy  was  chosen  solely  on  empirical  grounds.   Next* 
procedure  H  obtains  estimates  of  the  values  at  the  base  of  the  partial 
derivatives  of  each  violated  constraint  function  with  respect  to  each 
independent  variable.    Then  a  loop  is  initiated  which  attempts  to  obtain 
one  exploratory  move  for  each  presently-violated  constraint. 

It  is  now  necessary  to  return  to  the  subject  of  variable  coupling 
during  tangent  exploration.   The  manner  in  which  secondary  variables  are 
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selected  is  also  partly  shown  in  Chart  6.   It  should  be  recalled  that  the 
b-th  variable  is  decremented  by  an  approximation  of  the  amount, 


db  = 


^7 


D 


when  the  a-th  variable   is   primary   for  the   exploratory  move   in   a  hyperplane 
approximately  tangent   to  the  hyper-surface,    g-   =   0.      The   adjustment   of 

el 

I  by  d,  is  designed  to  confine  the  trial  move  to  that  hyperplane  when  X 

is  incremented  by  D  . 

a 

In  order  to  avoid  possible  computational  difficulties  arising  from 

large  disparities  in  magnitude  between  dgi  /ox„   and  og.  /O'x,  ,  it  was 

«j  a  j  d 

decided  simply  to  choose   a  and  b   such  that 


^g.  /^xb 


Y 


Here  the  notation 


refers  to  the  absolute  values  of  the  partial  derivatives.   Both  are 
evaluated  at  base  point  Y.    To  facilitate  this  rule  the  matrix  I  is  pro- 
duced by  auxiliary  procedure  H.   It  should  be  noted  that  the  partial 
derivatives,  stored  in  the  matrix  P,  are  correspondingly  arranged  at  exit 
from  H. 

Many  of  the  elements  of  P  will  be  zero  for  practical  problems.   In 
fact,  a  column  of  P  representing  an  explicit  constraint  (upper  or  lower 
bound  on  an  independent  variable)  contains,  at  every  violation,  (n  -  l) 
zeros,  while  the  last  element  is  always  unity.   Implicit  and  explicit 
constraints  are  not  treated  separately  by  the  algorithm  described  here. 
It  would  probably  be  advantageous  to  follow  the  example  of  M.  J.  Box  (3) 
and  do  so  in  future  work. 

When  some  element,  P  . ,  of  matrix  P  is  zero,  the  corresponding  i-th 

K.  J 

variable,  where  i  =  I   ,  cannot  be  used  either  as  a  primary  or  secondary 

kj 
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variable  in  the  sense  previously  defined.   In  that  case  the  tangent  hyper- 
plane  is  parallel  to  the  i-th  coordinate  axis.   Such  a  variable  will  be 
called  "passive"  in  the  following  discussion.   Any  variable,  the  q-th, 

where  q  =  I   ,  for  which  the  corresponding  P,  .  is  not  zero,  will  be  called 
kj  KJ 

"active". 

Since  any  increment  (or  decrement)  to  X.  cannot  remove  it  from  the 
tangent  hyperplane,  the  trial  step  for  any  passive  variable  can  be  chosen 
as  the  corresponding  element  of  the  D  vector.   Also,  if  more  than  one 
element  of  P,  . ,  (k  =  1,  2,  ...,  n)  are  zero,  any  two  such  passive  variables 
may  be  simultaneously  incremented  and/or  decremented  without  departure 
from  the  tangent  hyperplane.   In  that  sense  primary  and  secondary  passive 
variables  may  also  be  chosen  arbitrarily. 

Therefore,  in  the  selection  of  the  variables  for  the  trial  components 
of  the  tangent  exploratory  move  at  the  j-th  constraint,  the  primary  index 
can  always  become  in  turn, 

a  =  Xij'  Z2j,  .-•,   X(n  _  i)j' 

The  variable  of  index  I  •  is  never  primary.   The  corresponding  element 
PRi  cannot  be  zero,  since  this  would  imply  that  all  P^-i-  (k  =  1,  2,  ...,  n) 
are  zero,  in  which  case  the  tangent  hyperplane  would  be  parallel  to  all 
coordinate  axes.   Therefore,  it  cannot  be  involved  in  coupling  passive 
variables.   Also,  it  is  never  used  as  the  primary  in  coupling  active 
variables,  since  it  will  have  already  been  considered  at  least  one  time 
as  a  secondary  variable  by  the  time  a  =  1/     i ) • ' 

The  manner  in  which  the  secondary  variable  index,  b,  is  chosen  will 
be  termed  the  "coupling  mode".   At  all  times  the  j-th  element  of  N  con- 
tains the  current  coupling  mode  number  to  be  next  applied  in  calculating 
a  trial  move  in  the  tangent  hyperplane  at  the  j-th  constraint.   If  the 
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trial  exploration  results  in  a  tangent  move,  the  coupling  mode  is  unchanged. 
Otherwise  the  present  mode  number  is  decreased  by  one.   If  all  fail  in  turn, 
a  total  of  2(n-l)  modes  are  tried  to  complete  a  cycle.   If  necessary  it  is 
then  re-initiated. 

The  coupling  constants,  g  and  h,  for  each  constraint  are  functions  of 
the  current  corresponding  coupling  mode  number.   Constant  g  is  used  by 
procedure  T  to  couple  active  variables,  and  h  is  used  to  couple  passive 
variables.   Further  information  will  be  given  in  the  description  of  that 
auxiliary  procedure.   An  example  of  a  complete  cycle  of  g  and  h  is  given 
in  the  notes  to  Items  (8),  (9),  and  (10)  of  Chart  6. 

After  each  exit  from  procedure  T,  the  result  is  tested.   The  first 
trial  which  has  been  found  feasible  by  T  and  also  produces  a  value  of  the 
objective  function  below  the  base  value  becomes  a  "tangent  move".   The 
base  is  transferred  to  this  point,  and  the  vector  defined  by  that  transfer 
becomes  the  information  upon  which  a  pattern  move  is  attempted.   Exit  to 
that  part  of  the  main  procedure  immediately  follows. 

However,  if  after  trials  for  all  violations  are  complete  and  no 
tangent  move  is  made,  transfer  of  the  base  may  still  be  performed  by  a 
jump  move,  provided  certain  previously-described  conditions  are  met. 
Such  a  new  base  must  be  feasible,  even  though  it  may  involve  an  increase 
in  the  objective  function  over  that  of  the  previous  base.   This  transfer 
of  base  does  not  provide  a  promising  vector  for  a  pattern  move.   Instead, 
tangent  exploration  from  the  new  base  follows  any  jump  move. 

In  analyzing  Chart  6,  it  is  important  to  note  that  every  feasible 
trial  component  for  all  violated  constraints  is  tested  for  corresponding 
objective  function  value.   The  best  move  of  this  set  is  preserved  in  Z, 
the  function  value  in  f  ,  and  constraint  function  values  in  G.   Auxiliary 
procedure  T  indicates  that  at  least  one  feasible  trial  was  found  by  setting 
the  feasible  move  indicator,  t,  tcu unity. 
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Chart  6:   Preamble  to  a  Series  of  Tangent  Exploratory  Moves 
and  Tests  of  Results. 

Main  Procedure,  Part  III 
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Descriptive  titles  and  further  notes  for  the  items  of  Chart  6  are 
given  below. 

( 1 )  Is_  the  failure  counter  greater  than  zero? 

Here  the  counter  is  used  to  indicate  whether  the  preamble  has  been 
entered  after  failure  of  an  exploratory  or  a  pattern  move.   If  the  former, 
explore  step  sizes  are  reduced. 

( 2 )  Reduce  step  sizes ;  set  failure  counter  to  unity. 

( 3 )  Internal  auxiliary  procedure  H:v,  B,  G,  K,  X,  Y;  I,  P  is  performed. 
The  detailed  steps  of  this  auxiliary  are  shown  in  Charts  9  and  10 

below. 

( h )   Set  feasible  move  indicator  to  zero;  initialize  loop  to  be  indexed 
.1  =  1  (1)  v. 

( 5 )  Obtain  index  of  the  j-th  violated  constraint;  reset  trial  variables 
and  function  to  base  values. 

(6)  Is_  present  coupling  mode  number  zero? 
If  so,  cycle  must  be  restarted. 

(7)  Initiate  coupling  cycle. 

( 8 ) ,  ( 9 ) >  ( 10 )   Calculate  active  and  passive  coupling  constants. 
In  all  cases, 

g  =  N^  -  n  +  1  _or  Nj_ ,  whichever  is  positive, 
h  =  hQ  -  q ,  where  h   is  the  previous  value  of 
h,  and 

jql   =  n  -  g  -  1,  to  be  prefixed  by  the  sign  of  h  . 
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Therefore,  when  n  =  k,    a  complete  cycle  is 


Order 

N. 

1 

g 

h 

1st 

6 

3 

0 

2nd 

5 

2 

-1 

3rd 

h 

1 

+1 

Uth 

3 

3 

-2 

5th 

2 

2 

+2 

6th 

1 

1 

-3 

The  manner  in  which  g  and  h  are  used  to  select  secondary  variables 
is  described  in  the  discussion  of  auxiliary  procedure  T  below. 

(11)  Exploratory  Move  in  Tangent  Hyperplane . 

The  details  of  this  internal  auxiliary  procedure  are  shown  in 
Charts  7  and  8. 

( 12 )  Is_  trial  function  lower  than  value  at  base? 

If  so,  an  immediate  pattern  move  is  made  and,  in  this  process,  the 
base  is  transferred  to  the  tangent  move  just  found. 

(13)  Decrease  coupling  mode  number. 

4 
'  -Since  the  mode  used  by  procedure  T  at  the  i-th  constraint  did  not 

produce  a  tangent  move,  the  next  attempt  at  that  constraint  will  try  the 
next  coupling  mode. 

( lU )   Have  explorations  been  tried  in  all  hype rp lanes  tangent  to  presently- 
violated  boundaries? 

( 15 )  Increment  index  and  return  to  consider  next  constraint . 

(16)  Increment  failure  c ount e r . 

No  tangent  move  has  been  found.   Therefore,  return  to  the  main 
procedure  is  made  to  decrease  step  sizes,  test  for  end  of  search  and, 
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possibly,    a  jump  move. 

Special  Expedient   for  Convex  Constraint  Hyper-surfaces. 

It  should  be  noted  that   convex  curvature   of  a  violated  constraint  hyper- 
surface   is   inherently  troublesome   for  the   algorithm  described  here;    any 
move  which  were   strictly   confined  to   a  hyperplane   actually  tangent  to  a 
convex  constraint  hyper-surface  would  cause   immediate   re-violation  of  that 
constraint. 

Therefore,    a  special  expedient   is   employed  whenever  any  component   of 
a  trial  exploratory  move  produces   re-violation   of  the   constraint   currently 
considered.      This   open  sub-procedure   is   located  within  the   auxiliary  pro- 
cedure T.      This   problem  was   dealt  with  by  Glass   and  Cooper   (2)   by  their 
"extension"  to  the   alternate  move. 

The  method  proposed  here   takes   advantage   of  the   fact   that   a  measure 
of  the   gravity  of  a  subsequent   re-violation   is   provided  at  the   same   time 
it   is   detected.      When,   for  the   secondary  variable,    indexed  by  b, 

pb*Pgj    ^xi     Y     and 

G.<  0,  after  application  of  the  external  auxiliary  C(X),  the  following 
J 


^b  =  Xb   "  2Gj    /  Pb' 


simple   replacement   is   used  to  adjust  the  trial  move 

Xb  = 
A  factor  of  2   is   used  to  over-compensate   for  the   estimated  convexity.      If 

re-violation   follows,   this    adjustment   is   again  employed.      However,    if  a 
third  attempt   fails,   the  present   trial  component    is   abandoned.      Either  the 
next  primary  variable   is   selected  or  the   reflection   is  then  tried. 

The   steps   outlined  above   are   appropriate   only  when  the  variables    in- 
volved are   active.      If  the   primary  variable   is   passive   and  not  presently 
coupled,    the   variable   of  largest  M^g.    /^x.[         I    ,    (i   =  1,   2,    ...,   n), 
is   selected  to  serve   as   the   secondary  variable   for  this   special  expedient. 
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The   index  of  that   variable   is   the   element   I    . . 

Finally,    if  two  passive   variables   are   coupled  and  re-violation   of  the 
currently-considered  constraint    occurs    on   any  trial   component,   the   second- 
ary variable   of  that    component  becomes   the  k-th  where 

Then   all   calculations    for  the    component    are   restarted.      Detailed  steps 
of  this   expedient   are    shown   in   Charts   7   and  8  below. 

Tangent  Exploratory  Auxiliary  Procedure. 

Input  parameters   of  this   procedure,   T,    include   g  and  h,   the   coupling 
mode   constants   of  the   current   trial  move;    i    and  j,    indices   of  the    con- 
straint hyper-surface  to  be   considered;    the   column   I    .,    (k  =  1,  2,    ...,   n), 
which  holds   indices   of  variables   ordered  by  magnitude   of  the    corresponding 
P      ,    (k  =   1,    ...,   n);    and  the   column  matrix  G,    containing  the   constraint 
function  values   at  X,  where   one   or  more   constraints   are   violated. 

This   internal  auxiliary  procedure  usually   generates  values    for  t,    f   , 
A,    and  Z,    all   of  which  are   utilized  only  if  no  tangent  move   is    found  for 
any   one   of  the    currently-violated  constraints.      Values   of  f  and  X,   pro- 
duced by  T,    constitute   a  tangent   move   if  and  only   if   f  lies   below   fy,  the 
value   of  the   objective   function   at   the  base  point,   Y. 

An   introduction  to  the   steps   of  this   procedure  has  been   given   in  various 
previous   sections.      However,    it   is   still  necessary  to   state  the   rules  by 
which  g  and  h  are  used  to  select   secondary  variables.      If  g  is  the   coupling 
mode   constant   for   active   variables   and   a  =   Iv-p   then  b   is    chosen   as   the 
(k  +  g)th   index   of  the   j-th   column   of  I,   except   that  the   index  at   In* 
pre-empts   in   all   cases  when  k  +   g  >  n. 

Similarly,    if  h   is   the   coupling  mode   constant    for  passive   variables 
and  a  =   Ijh  ,   then  b   is   chosen   as   the  w-th   index   of  the   j-th   column 
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whenever  w  =  k  +  I  hi  £  n-1.   If  w  >  n-1  or  if  P^-/  0,  no  coupling  is 
effected.   When  the  increment /decrement  of  the  b-th  variable  is  actually 
used,  the  sign  of  this  operations  is  dictated  by  the  sign  of  h. 

As  an  example,  we  shall  assume  that  n  =  U,  m=  3,  v=2,  K-,  =  2,  Kg  =  3, 
and 


I  = 


3  1 

1  2 
1+  k 

2  3 


P  = 


.1 

0 

-.5 

0 

7.1 

0 

-9.8 

1 

Then  the  variables  involved  in  the  trial  components  of  the  tangent 
exploratory  move  near  the  hyper-surf ace,  g2  (X)  =  0,  would  be,  for  the  first 
three  coupling  modes: 


N. 

l 

g 

a 

b 

6 

3 

3 

2 

1 

2 

h 

2 

5 

2 

3 

k 

1 

2 

1+ 

2 

k 

1 

3 

1 

1 

h 

h 

2 

The  last  three  modes  of  the  active  coupling  cycle  repeat  the  above  table. 
Since  no  element  of  the  first  column  of  P  is  zero,  only  active  coupling 
is  used  at  this  constraint.   The  variable  coupling  used  in  trial  components 
near  g^(X)  =  0  is  shown  in  the  table  below.   This  constraint  boundary  is 
actually  a  hyperplane,  perpendicular  to  the  3rd  coordinate  axis.   Therefore, 
all  coupling  is  passive. 
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N. 

h 

a 

b 

Remark 

6 

0 

No 

coupling 

a  =  1,  2,  U 

i 

5 

-1 

1 
2 
1+ 

2 

None 

Increments /Decrements 
opposite  in  sign 

k 

+1 

1 
2 

2 

u 

None 

Increments /Decrements 
the  same  in  sign 

3 

-2 

1 
2 

None 
None 

Increments /Decrements 
opposite  in  sign 

2 

+2 

1 
2 

None 
None 

Inc  rement  s /De  creme  nt  s 
same  in  sign 

1 

-3 

No 

coupling 

a  =  1,  2,  k 

Detailed  steps  of  procedure  T  are  shown  in  Charts  7  and  8.   De- 
scriptive titles  for  the  items  given  and  additional  explanatory  notes 
follow  the  charts. 
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Chart  7:   Internal  Auxiliary  Procedure  T:  g,h,i,j,G,I;  f,f  ,t,A,X,Z. 
The  Tangent  Exploratory  Move,  Part  I 
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Chart   8:      Internal  Auxiliary  Procedure  T:    g,h,i,j,G,I;    f,f    ,t,A,X,Z. 
The  Tangent   Exploratory  Move,   Part   II 
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Charts   7   and  8 

(1)  Initialize   index   of  outer  loop,   k  =   l(l)n-l. 

Index  k  will  pick  up,    in   order,   the   variable  numbers   of  the  n-1 
primary  variables   of  the  present   tangent   exploratory  move. 

( 2 )  Select   index   of  primary  variable ;    select  the   corresponding  partial 

derivative;    save   old  value   of  primary;    increment   and  save  new 

primary  variable  value . 

The   values   of  p    ,   y„ ,    and  x      are  preserved  separately   from  the 
a        °*  a 

matrices   in  which  they  are   found.      This  policy  reduces    indexing  steps    in 
an   object  program.      It    also   simplifies   the   calculation   of  the   component 
reflection,    if  this  becomes   necessary. 

(3)  Is   the  primary  partial  derivative,     jg.    /  fix      zero? 

If  so,    any   change   in  the   a-th  variable  would  be  parallel  to  the 
tangent  hyperplane   and,   therefore,    it   is  passive, 
{h)      Select   the   index,   b_,    for  the   secondary  active   variable  by  using 

constant  g_;    select   secondary  partial  derivative;    save   old  value 

of  the   active   secondary;    save  new  decremented  variable  value . 

Values   of  p-,  ,   y,  ,    and  x-.     are   also  preserved   for  later  use   in 
calculating  a  trial  reflection. 
( 5  )      Change  effective   value   of  secondary   variable  ,   X-u, . 

(6)  Does   coupling   constant   indicate  that   the  primary  passive   variable 
is   to  be   coupled? 

If  so,    it  will  be  necessary  to   select   the   appropriate    index   of  a 
passive   variable. 

( 7 )  Compute  tentative   location   of  the   secondary   index . 

(8)  Is_  location  within  bounds? 

If  so,    It   is   still  necessary  to  discover  whether  or  not   the  tenta- 
tively-selected variable   is   also  passive. 


( 9 )      Save  the  partial  derivative  of  tentative   secondary  passive  variable. 

( 10 )  I_s_  the  partial  derivative   zero? 

If  not,    coupling  is  not  allowed  since   any  change   in  one  active 
variable  would  cause   departure   from  hyperplane. 

(11)  Select   secondary  passive   index,   b_,    from  I_  matrix;    save   old  value   of 
X,  ;    compute   incremented  or  decremented  secondary  variable. 

The   direction   of  that   adjustment   depends   upon  the  present   sign  of  the 
coupling  constant. 

( 12 )  Select  variable  with   largest   absolute  partial  derivative   as   secondary; 
save   value   of  derivative   and  old  value   of  variable;   maintain  this 
value   for  trial. 

These   steps   allow  the   reflection/reset   logic   to  handle  non-coupled 
components  without   special  tests.      They  also  provide   a  secondary  variable 
for  possible   use   in  the   special  expedient   for   convex   surfaces. 

(13)  Set   re  fie  c t i  on /re  set   switch  to  unity. 

This   switch   is  tested  whenever  a  trial  component   is  non-feasible 
or   fails   to   improve  the   objective   function  value.      When  s   is   unity,   the 
reflection   of  the  present   trial  component   is   next   tested;    if  s   =  2,   the 
values   of  the   variables  which  existed  prior  to  the   failed  component    are 
reset. 
( 1^ )      Initialize   loop   to   count   possible    applications   of  the   special 

expedient   at   a  convex  hyper-surface . 

This   sub-procedure   is   tried  a  maximum  of  three   times   consecutively 
before   a  trial  component   is   abandoned. 
( 15 )      Compute   constraint    functions   and  test    for  violations. 
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(16)  Is_  any  constraint  violated? 

If  not,   exit   from  the   special  expedient   loop   is  made.      In  that  case 
the  next   step   is  to  test  whether  or  not   this    feasible   trial   component 
brings   an   improvement   in  the   objective   function  value. 

( 17 )  Is_  the  presently-considered  constraint   re-violated  by  this  trial 
component? 

If  not,   the   special  expedient   is   useless   and  this   trial  component 
has   failed. 

(18)  Is_  the   secondary  variable  passive? 

If  so,  it  cannot  be  used  in  the  logic  of  the  special  expedient. 
Therefore,  exit  is  made  in  order  to  try  the  reflection  or  to  select  a 
different  secondary  variable. 

If  not,   the   secondary  variable   is   adjusted  by  approximately  twice 
the  minimum  amount  necessary  to  produce   a  feasible  trial.      These   steps 
are   discussed  at   Item  (26),   etc.,  below. 

( 19 )  Calculate  the   value   of  the   objective   function   corresponding  to 
the  present   trial  component . 

(20)  and   (21)      Is   this  value   an   improvement   over  the  best  previous 

result    at  this    constraint? 

If  so,    it    is   saved,   the  variables   remain  undisturbed, 
and  the  next  primary  variable   is    selected. 

( 22 )  Has   any   feasible  move  been  previously   found  at    any  presently- 
violated  constraint? 

If  not,   this    component  will  be  preserved. 

(23)  Set   feasible  move   indicator  to  unity. 

(2k)      Save  present   values   of  objective   function,    independent  variables, 
and  constraint   functions   as  best   feasible  move   so  far. 
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(25)  Test  objective  function. 

The  function  value  at  the  present  trial  component  is  compared  with 
the  value  for  the  best  feasible  component  so  far  at  all  constraint  surfaces 
currently  violated. 

If  it  is  below  the  previous  best  value  the  appropriate  quantities 
are  preserved  for  a  possible  jump  move.   If  not,  the  reflection/reset 
switch  is  tested  to  determine  the  manner  of  finding  the  next  trial  com- 
ponent.  See  Item  (29),  below. 

(26)  Correct  secondary  variable. 

(27)  and  (28)   Have  three  attempts  to  apply  the  special  expedient  been 

made? 
(29)  and  (30)   Test  reflection/reset  switch. 

If  a  trial  component  fails,  its  reflection  in  the  hyper- 
plane  about  the  previous  best  reference  is  tried.   If  a  reflection  fails, 
variables  are  reset  to  the  previous  reference. 

The  switch  test  at  Item  (30)  is  made  after  failure  of  the 
special  expedient. 

( 31 )  Compute  reflection  of  trial  component  in  hype rp lane  about  previous 
best  point;  set  re  fie  cti  on  /res  et  switch  to  2_. 

It  should  be  noted  that  the  reflection  of  a  trial  which  is  modified 
by  the  special  expedient  is  not  used.   Control  next  returns  to  the 
feasibility  test. 

( 32 )  Reset  primary  and  secondary  variables  to  previous  best  values . 
Both  the  original  and  reflection  components  for  this  pair  of 

variables  have  failed. 

(33)  and  (3*0   Have  (n-l)  primary  variables  been  tried? 

If  so,  this  application  of  T  is  complete.   If  not,  k  is 
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incremented  and  the  next  primary  variable  is  selected. 
( 35 )   Reset  secondary  variable. 

If  both  components  fail  with  passive  coupling  due  to  convexity  of 
the  constraint  hyper-surface,  a  new  secondary  variable,  which  is  neces- 
sarily active,  is  selected. 

Calculation  of  Partial  Derivatives  of  the  Constraint  Functions. 

In  the  material  above  discussion  of  the  problem  of  evaluating  the 
partial  derivatives, 


,   (J  =  1,  2,  ...,  v),  (i  =  1,  2,  ...,  n), 

Y 


has  been  bypassed.   In  Chart  6  this  process  is  shown  as  an  internal 
auxiliary  procedure 

H:  v,  B,  G,  K,  X,  Y;  I,  P   , 
which  is  used  as  part  of  the  preamble  to  a  series  of  tangent  exploratory 
moves.   The  details  of  one  formulation  of  that  auxiliary  are  shown  in 
Charts  9  and  10.   This  is  the  algorithm  used  in  the  test  calculations 
to  be  discussed  in  the  following  sections. 

However,  there  are  many  computational  alternatives  for  accomplishing 
the  same  purposes,  one  or  more  of  which  may  be  more  elegant  and  economical 
than  that  given  in  Charts  9  and  10.   The  last  section  of  this  paper  will 
return  to  this  point. 

Five  major  presuppositions  led  to  the  particular  version  of  procedure 
H  which  is  described  here. 

1.   Only  one  external  auxiliary  procedure  to  define  all  constraint 
functions  should  be  allowed.   All  such  functions  should  be  evaluated  each 
time  this  auxiliary,  C(X),  is  called. 

The  author  has  visualized  C(X)  as  a  subroutine  to  be  supplied  by 
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the  user  of  a  standard  library  routine  which  contains  all  steps  except 
F(X)  and  C(X).   For  this  reason  complete  simplicity  of  logic,  e.g.,  no 
conditional  branching,  within  both  of  these  external  auxiliaries  has  been 
an  overriding  desideratum. 

2.  Wo  special  external  auxiliaries  to  calculate  the  derivatives, 
themselves,  either  by  difference  methods  or  by  the  evaluation  of  formulas 
obtained  by  analytic  differentiation,  should  be  allowed. 

The  present  author  felt  that  the  entire  task  of  obtaining  de- 
rivatives should  be  shouldered  by  any  proposed  minimization  method.   This 
is  particularly  important  for  any  method  seriously  suggested  for  the 
solution  of  practical  problems. 

3.  Special  software  or  subroutines  to  obtain  derivatives  by  the 
methods  of  Wengert  (5)  or  Smith  (6)  should  not  be  assumed. 

k.      Implicit  and  explicit  constraints  would  not  be  separately  treated. 

5.   The  minimum  number  of  evaluations  of  F(X)  and,  expecially,  C(X) 
in  obtaining  satisfactory  solutions  would  be  a  primary  goal. 

In  the  context  of  these  presuppositions  it  was  decided  to  approximate 
the  desired  partial  derivatives  by  simple  first-order  differences.   In 
other  words  the  secant  is  actually  substituted  for  the  tangent. 

When  procedure  H  is  called, one  or  more  tangent  exploratory  moves  are 
about  to  be  tried.   The  actual  number  is,  at  most,  v,  the  number  of  con- 
straints violated  by  the  last  trial  move.   The  point  Y  is  the  present 
base,  to  which  corresponds  the  column  matrix  of  constraint  function  values, 
B.   All  elements  of  B  are  zero  or  positive.   The  point  X  is  the  non- 
feasible  trial  move,  to  which  corresponds  the  column  matrix,  G,  v  elements 
of  which  are  negative.   Additionally,  K  is  a  column  matrix  of  v  elements, 
each  of  which  is  the  index  number  of  a  violated  constraint. 
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In  the  following  discussion  attention  will  be  restricted  to  the 
approximation  of 

,   (1-1,  2 n), 


^Xi 

where  k  is  any  one  of  the  indices  stored  in  K,  and  g  is  the  k-th  constraint 

k 

function.   It  should  be  understood  that,  each  time  it  is  called,  H  produces 
all  elements  of  the  matrix,  P,  which  contains  n  rows  and  v  columns. 
The  total  differential  of  e  can  be  expressed, 


where  the  x.   (i  =  1,  2,  ...,  n)  are  general  independent  variables. 
Evaluation  of  the  partial  derivatives  is  desired  either  at  the  base  point 
or  nearby,  for  example,  at  some  point  between  the  current  base,  Y,  and  the 
non-feasible  trial,  X.   The  differentials  above  are  approximated  by  simple 
differences , 

dx.  ^  X.  -  I.,     (i  =  1,  2,  ...,  n), 
and 

d^~  Gk"  V 

It  is,  in  general,  necessary  to  evaluate  only  (n-l)  of  the  partial  de- 
rivatives ^gv/^x.,  (i  =  1,  2,  ...,  n),  since  the  n-th  can  be  found  from, 


L 


ZllVi 

i=l 
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One  call  of  the  external  auxiliary  procedure  C(X)  is  necessary  to 
calculate  each  d  gk/9x.  for  *■  =  l(l)n-l.   These  values  are  then  approxi- 


1 
mated  by  the  simple  differences, 


^gk/3Xi2:  (Ak  -  B  )  /  (X.  -  Y±), 


where  A,  is  g,  (W) . 

The  point  W  is  defined, 

W  =  Y  ,  for  all  q^  i, 
q    q 

and  W-?  =  X.  . 

±  1 

The  method  described  above  would  fail  whenever  any  X.  =  Y.   When  • 

this  happens  W^  is  taken  as 

W.  =  X.  +  D.  , 
1    l    l' 

in  the  evaluation  of  A^. 

In  the  notation  of  the  flow  charts 


^  P..  . 

y        lk 


"}  V^Xi| 

Auxiliary  procedure  H  generates  matrix  P  containing  n  rows  and  v  columns 

of  such  approximate  partial  derivatives  corresponding  to  the  n  variables 
and  v  violated  constraints.   At  exit  from  H  each  column  of  P  is  arranged 
in  ascending  order  of  absolute  value.   The  procedure  also  generates  matrix 
I  which  contains  the  variable  numbers  corresponding  to  the  arrangement  of 
P.   The  sorting  which  is  necessary  to  order  P  and  I  is  accomplished  by  the 
simple  exchange  method. 

The  steps  which  calculate  the  elements  of  matrix  P  are  given  in  Chart  9- 
The  rearrangement  of  P  and  generation  of  I  are  shown  in  Chart  10.   The 
author  has  omitted  the  usual  listing  of  descriptive  titles  for  these  charts, 
since  there  are  no  difficult  points  requiring  more  detailed  explanation. 
It  should  be  mentioned  that  the  control,  e,  and  switch  s  are  used  for 

hi 


branching  when   any  increment,    Z.    =    (X-    -  Y.  )  ,    is   near   zero   and  t< 
eliminate   the   unnecessary  use   of  C(X),    as   explained  above. 
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k=K. 


Pij    "   <VBk)/x 
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Chart  9:  Internal  Auxiliary  Procedure  H:  v,B,G,K,X,Y;  I, P. 

Evaluation  of  Partial  Derivatives  of  Constraint  Functions, 
Part  I 
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© 


Ii) 


p  .  -  p  Vz 
qj   qj  q 


(3) 


(1-1,2,. ...n) 


0 
I 


JL 


\j      (4) 


19) 


>  k  -  k+1 


-KD 


(11) 


Chart  10:   Internal  Auxiliary  Procedure  H:  v,B,G,K,X,Y;  I, P. 

Evaluation  of  Partial  Derivatives  of  Constraint  Functions, 
Part  II 
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Numerical  Results  and  Comparisons 

The  new  calculations  reported  here  were  performed  on  the  Control  Data 
Corporation  l60U  Computer  of  the  Computer  Facility,  U.  S.  Naval  Postgraduate 
School,  Monterey.   Both  the  Tangent  Search  Method  described  above  and  the 
Complex  Method  of  M.  J.  Box  (3)  were  programmed  by  the  author  in  FORTRAN  60 
for  these  trials. 

Comparisons  between  other  methods  and  Tangent  Search  are  made  on  the 
basis  of  criteria  used  by  previous  authors.   Counts  of  required  moves  are 
given  for  comparison  with  the  results  of  Glass  and  Cooper  (2),  and  for  the 
one  problem  reported  in  detail  by  Klingman  and  Himmelblau  (l).   However, 
most  comparisons  are  between  Tangent  Search  and  the  Complex  Method.   All 
eleven  examples  were  run  with  both  of  these  methods  by  the  author.   These 
comparisons  are  made  in  terms  of  the  numbers  of  function  and  constraint 
evaluations,  tabulated  separately,  required  to  attain  a  constrained  minimum 
of  a  defined  accuracy.   The  present  author  wants  to  concur  in  the  remark 
of  M.  J.  Box  (3)  concerning  the  superiority  of  this  measure.   However  for 
Tangent  Searches  counts  of  base  points,  conventional  exploratory  moves, 
and  tangent  exploratory  moves  are  also  tabulated  below  for  possible  help 
to  other  investigators.   For  the  Complex  Method,  total  and  permissible 
trials  are  counted  in  addition  to  evaluations  of  the  external  auxiliary 
procedures.   These  statistics  are  used  by  M.  J.  Box.   It  is  the  interpret- 
ation of  the  present  author  that  a  trial  vertex  is  "permissible"  only 
if  it  both  is  feasible  and  also  produces  a  value  of  the  objective  function 
which  is  no  longer  worst. 

Each  trial  problem  will  be  defined  below  and  results  tabulated.   A 
summary  of  the  results  is  given  subsequently. 
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PROBLEM  1 

This  example  is  taken  from  the  paper  of  Glass  and  Cooper  (2).   Their 
Table  1  shows  the  solution  as  developed  by  Sequential  Search. 

Required  is  the  minimum  of 

F   (x)  =  -  V  21  -  (Xi  -  5)2  -  (x2  -  5)2, 

subject  to  the  constraints, 


and 


G1  (x)  =  x±2  -   h   x2  2  0, 

G2  (x)  =  (x2  -  6)2   -  h(x±   -  3)  21  0  . 
All  trial  runs  were  calculated  from  the  initial  point, 

x  =  7.000 

x2  =  1.000 

F  =  -2.236,   with 
initial  step  sizes, 

D1  =    .0666667  and 

D2  =    .1333333  . 

According  to  the  article  the  actual  solution  to  the  problem  is 

x_   =  h. 

Is 

c2s 

F  =  -U. 7958  . . . 

s 

Sequential  Search  attained  a  result  within  a  circle  of  radius  .0015 
about  the  true  solution  as  its  23rd  base  point.   Tangent  Search  solutions 
were  also  obtained  to  the  same  accuracy  for  three  choices  of  step  reduction 
factor,  r.   Comparison  of  all  these  results  is  given  in  Table  la. 
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TABLE   la 


End  criterion: 


~]j(x1  -  kf  +   (x2  -   k)2        6     .0015 


Result 

Algorithm 

Seq.    Search 

Tangent   Search 

r  = 

"     -25 

.375 

.500 

xl 

I2 

F 

U.000 

3.999 

-U.796 

U.0005 

3.999^ 

-^•7958 

3.999^ 
3.9988 

-^•7955 

3.9992 
U.0008 

-^.7958 

Base 
Points 

23 

21 

31 

32 

Explore 
Moves 

12 

5 

7 

11 

Alternate 
Moves* 

5 

6 

12 

19 

*     For  the  purposes   of  Tangent   Search  this   counter  was    incremented  before 
any  entry  to  Item  5   of  Chart  6. 

Results  with  Sequential   or  Tangent  Search   seem  substantially  the   same 
for  this   problem.      The  best   run   is    for  Tangent  Search  with  r  =    .25    . 
However  with   larger   factors   this  method  attains   the   required  minimum  less 
efficiently  than  Sequential  Search. 

Problem  1  was    also   solved  using  the   Complex  Method  of  M.    J.    Box    (3)    • 
Initial   conditions  were,    as  before,   x-.    =   7.0   and  x     =   1.0    .      Two  runs  were 
made   in   order  to  discover,   to  a  slight  extent,   the  effect   of  the   random 
choices   of  vertices  made   at   the  beginning  of  each   calculation.      The   sequence 
of  pseudo-random  numbers,   used  in   setting  up  the   first   complex,   varies 
between  the   two  solutions.      These   are   denoted  "a"   and   "b".      To   compare 
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performance  with  Tangent  Search  the  number  of  functions  and  constraint 
evaluations  were  counted.   For  the  Complex  Method  the  centroid  of  the 
vertices  is  the  reference  point  for  determining  the  progress  of  the  solution. 
Comparison  is  made  for  the  same  degree  of  accuracy  described  previously; 
results  are  given  in  Table  lb. 

TABLE  lb 


End  Criterion: 


V  (^  -  M2  +  (x~  -  h)2     i-    . 


0015 


Algorithm 

Number  of  Evaluations 

Complex  Method 

Constraints 

Ob  j .    Function 

a 
b 

191 
128 

109 
61+ 

Tangent   Search 

k2 

38 

r  =    .125 

r  =   .375 

63 

72 

r  =   .500 

81     . 

85 

These   results   seem  to  indicate  that  Tangent   Search   is   superior  to  the 
Complex  Method  for  this  problem. 
PROBLEM  2 

This   is   the   second  illustration  of  Glass    and  Cooper.      The  problem  is 
to  minimize   the   objective   function   of  Problem  1,    above,   but   now  subject   to 
the   single   constraint, 

G(x)   =   32-Ux     -  x2a  >  0. 
The   initial  point   and  step   sizes    are,    as  before, 


x1  =   7.0    . 
x,   =   1.0   > 


D     =    .0666667, 
D2  =    .1333333. 


5U 


The  solution  (to  four  decimal  places)  is  said  to  be: 


xls  =  U.37U1 


xQ  =  3.8083 

2S 


F  =  -U.815U  • 

s 

Table  2  of  that  reference  shows  the  steps  generated  by  Sequential 
Search  until  a  base  approximately  .  00U  units  from  the  true  result  is 
found.   Previous  base  points  were  actually  closer  than  this.   However 
none  of  these  was  chosen  as  final  since  subsequent  results  diverged. 
Comparisons  between  that  result  and  corresponding  tests  using  Tangential 
Search  are  given  in  Table  I la. 

TABLE  Ila 


End  Criterion:       W  (x     -   x      )      + 


(x     -   x2s)2       6     .00k 


Result 

Algorithm 

Seq.    Search 

Tangent   Search 

r  = 

.25 

.375 

.500 

Xl 
x2 

F 

U.376U 

3.8071 

-U.815U 

U.37H1 

3.8081 

-U.815U 

U.37U1 

3.8083 

-U.815U 

U.37^1 
3.8082 

-U.815U 

Base 
Points 

hi 

36 

32 

36 

Explore 
Moves 

20 

10 

13 

15 

Alternate 
Moves 

15 

lit 

15 

21 

Since   two  runs    confirm  the   conclusion,    it   seems  that  Tangent    Search 
is  better  than  Sequential  Search   for  Problem  2.      Either  of  the   two  smaller 
values   of  r  seems   to  be   the  proper  choice   in  this   case. 
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Trial  calculations  were  also  made,  as  on  Problem  1,  with  the  Complex 
Method.   These  results  are  shown  in  Table  lib. 


TABLE  lib 


End 


Criterion:   ~W  (x1   -  *ls)2  +  U2  "  X2S^  -  -00^ 


Algorithm 

Number 

of  Evaluations 

Complex  Method 

Constraints 

Ob j .  Function 

a 

132 

68 

b 

122 

65 

Tangent  Search 

lh 

81 

r  =  .250 

r  =  .375 

72 

8U 

r  =  .500 

91 

101 

The  author  feels  this  table  shows  substantially  similar  performance 
for  the  two  methods . 

PROBLEM  3 

This  example  is  taken  from  the  paper  of  Klingman  and  Himmelblau  (l); 
however,  it  originated  in  another  article  by  R.  A.  Mugele  where  it  demon- 
strated the  IBM  Probe  Method  of  Optimization.   A  detailed  narrative  of  the 
progress  of  the  solution  by  the  Multiple-Gradient  Summation  Technique  is 
provided  by  the  Appendix  to  the  former  paper.   It  is  required  to  locate  the 
optimum  of 

F(x)  =  1/   ((x,  +  l)2  +  x  2), 


1 


under  the  two  constraints, 


G±   (x)  =  x  2  +  x22  -  k   Z    0 


G2  (x)  =  16  -  xx2  -  x22  >  0. 
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To  prevent  possible  confusion,  this  problem  will  be  immediately  re-cast 
to  require  the  minimum  of 

F(x)  =  -1/  ((x  +  D2  +  x22), 
subject  to  the  given  constraints. 

Initial  conditions  for  trial  runs  are, 

x1  =  5,  x  =  k9    and 

D1  =  1.5,  D2  =  1.2  . 

Since  G  (x)  =  -25  for  this  origin,  the  solution  must  begin  by  search- 
ing for  some  feasible  point.   Both  the  Multiple  Gradient  Summation  and 
Tangent  Search  Methods  include  a  preliminary  procedure  for  this  purpose. 
A  comparison  of  results  is  provided  by  Table  Ilia.   Moves  required  to  find 
an  initial  feasible  point  and  subsequent  performance  criteria  are  tabulated. 
In  order  accurately  to  compare  results  the  constrained  minimization  phases 
of  all  three  Tangent  Search  Calculations  were  started  at  the  point  (2.8, 
2.0)  at  which  all  counts  were  reset  to  zero.   Final  tabulations  are  made 
for  the  point  in  each  calculation  which  corresponds  in  accuracy  to  the 
last  point  reported  in  the  Appendix  of  (l).   Specifically  this  is  the  point 
at  which  the  base  enters  a  circle  of  radius  less  than  .52  about  the  true 
solution, 


Xls    = 

-2.0 

xls    = 

0 

F 
s 

-1.0 

Evaluation  of  analytic  first  partial  derivatives  are  also  counted 
for  Multiple-Gradient  Summation.   Other  criteria  are  as  described  pre- 
viously. 
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TABLE   Ilia 


End  Criterion: 


V^T 


+  2)2   +  x02    ±   .52 


Criterion 

^_ , 

Algorithm 

Initial 
Feasible   Pt. 

M-G  Summation 

Tangent   Search 

*•'*:= 

.250 

.375 

.500 

xl 

2.8 

1.625 

3.500 

-2.125 

Xo 

2.0 

1.300 

1.150 

-1.700 

F" 

-.059^ 

-.1165 

-.01+61+ 

-.2U06 

Base 
Points 

1 

2 

3 

3 

Explore 
Moves 

1 

1 

2 

1 

Constraint 
Evaluations 

6 

T 

15 

11+ 

Final  Point 

-2.1+9 

-2.118 

-2.101+ 

-2.131 

xl 

*2 

F^ 

-.182 

.1+699 

-.331+7 

.1+960 

-.1+005 

-.6799 

-.7513 

-.655^ 

M 

oves 

7 

15 

8 

8 

Base 

Explore 

10 

5 

1+ 

6 

Alt. 

3 

5 

1+ 

6 

Ev 

aluations 

1+3 

31 

20 

28 

Constraints 

Dbj.    Funct. 

33 

U3 

29 

1+3 

Part.    Deriv 

6 

— 



Inspection  of  Table  Ilia  does  not  seem  to  establish  any  superiority 
between  the  methods.   Among  step  reduction  factors  for  Tangent  Search, 
r  =  .375  seems  best.   With  that  choice  the  final  result  is  better 
than  that  of  Multiple-Gradient  Summation.   That  method  seems  to  do 
better  in  the  initial  search  for  a  feasible  origin;  however  the  author 
feels  this  is  due  to  a  particularly  fortuitous  initial  pattern  move. 
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Comparisons  between  Tangent  Search  and  the  Complex  Method  are  given 
in  Table  IITb.   The  end  criterion  is  smaller  for  these  calculations. 

TABLE  I lib 


End  Criterion: 


Y^T 


+2)      +  x^        fi    .005 


Algorithm 

Evaluations 

Trials/Moves 

Complex  Meth. 

Constraints 

Ob.i .  Funct . 

Total 

Permissible 

a 

295 

223 

287 

180 

b 

- 

- 

failed 

- 

T 

angent  Search 

Base 

Explore 

Alt. 

r  =  .250 

68 

87 

33 

11 

13 

r  =  .375 

71 

86 

30 

10 

16 

r  =  .500 

75 

101 

25 

15 

19 

Using  the   sequence   of  pseudo-random  numbers   designated  "b" ,   the   Complex 
algorithm  failed  to  develop  a  feasible   fourth  vertex  at   the   start.      This 
was   due  to  the   fact   that  both  the   random  trial  vertex  and  the   centroid  of 
the  three  previously-chosen  vertices  were  not    feasible.      Both   lay  within 
the   inner  boundary. 

G1(x)    =   x±2   +   xg2   -   1+   =   0. 

This  difficulty  may  indicate  an  undesirable  restriction  on  the  use  of 
that  method  in  some  practical  problems.   The  author  feels  that  this  failure, 
with  the  poor  result  for  run  "a",  indicates  the  superiority  of  a  pattern 
search  for  Problem  3.   Finally,  the  results  given  for  Tangent  Search  in 
Table  I lib  are  better  for  the  two  smaller  step  reduction  factors  considered. 
PROBLEM  k 

This  example  is  based  upon  the  "banana-shaped  valley"  proposed  origin- 
ally by  Rosenbrock  (7).   A  modification  of  this  function  was  used  in  the 

calculations  reported  by  Krolak  and  Cooper  (8).   Problem  k   calls  for  minim-. 

...  < 

lzmg  this   function  under  the   constraints   proposed  for  the   second  problem 
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of  Table  1  in  the  article  of  Klingman  and  Himmelblau  (l). 
Specifically,  the  minimum  of 

F(x,y)  =  (y  -  x2)2  +  (l  -  x)2 
is  required,  assuming  the  constraints, 

G1   =  x  -  .2  2    0 

G2  =  2  -  x  Z    0 

G3  =  y  -.2  Z.    0 

Gh   =  2  -  y  Z    0 

G5   =   1  -  x2  -  y2   Z     0. 
It  will  not  be  possible   to   compare  trial   calculation  with  the  Multiple- 
Gradient   Summation  Technique   since   it    failed  on  this  problem.      From  evidence 
in  the  text   it   is   assumed  that   Problem  U,   here   defined,    is   the   same   as   the 
second  problem  of  Klingman   and  Himmelblau,    and  that  the   function   actually 
shown   in  their  table   contains   a  typographical  error.      Exact   comparisons 
are   impossible   for  another  reason;    the   starting  points   of  all  of  the  tests 
reported  in  their  Table   1  are   omitted. 

Tangent   Search  Calculations  were   at   first   started  with  the  non-feasible 
point,    (5 s   -5)    and  with  both   step   sizes  unity.      Table   IVa  displays   totals 
of  moves   and  evaluations   required  to   find  feasible   origins,  which   are   also 

given. 

TABLE   I Vaj  Tangent   Search 


Reduction 
Factor 

Feasible   Origin 

Moves 

Constraint 
Evaluations 

X 

y 

F 

Base 

Expl. 

.250 
.375 
.500 

.333 
•  73U 
.250 

.766 
.kk-J 

.875 

.Qlh 

.079 

1.223 

6 
6 
h 

3 
5 
3 

17 
30 
13 
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In  order  validly  to  compare  step  reduction  factors  in  Tangent  Search 
and  that  algorithm  with  the  Complex  Method,  all  of  the  results  in  Table  IVb 
were  obtained  for  calculations  initiated  at  the  feasible  point,  (  .25,  .875) 
Beginning  step  sizes  for  Tangent  Search  runs  were  unity  in  all  cases.   The 
correct  solution  to  Problem  h   is  (to  six  significant  figures  in  each  inde- 
pendent variable ) , 

xs  =   .808169 

ys  =   .588951 

Fs  =   .0U09190. 


TABLE  IVb 


End  Criteria:  I  be  —  x 


k< 


000005    and 


f-4r 


000005 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Const . 

Ob.i .  Funct. 

Totals 

Permissible 

a 

179 

12  h 

169 

90 

b 

196 

131 

186 

96 

Tangent  Search 

61+ 

80 

Base 

Explore 

Alt. 

r  =  .250 

27 

Ik 

15 

r  =  .375 

96 

109 

37 

18 

2U 

r  =  .500 

98 

118 

33 

21 

26 

All  three  Tangent  Search  results  are  superior  to  both  Complex  Method 
solutions.   Additionally,  a  reduction  factor  of  .25   is  best  for  this  prob- 
lem. 
PROBLEM  5 

This  trial  problem  is  also  adapted  from  Krolak  and  Cooper  (8)  with 
constraints  proposed  by  Klingman  and  Himmelblau  (l).   It  is  the  minimi- 
zation problem  corresponding  to  the  eighth  example  of  their  table  1. 
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Specifically,  it  is  desired  to  find  the  minimum  of 

F(x)   =  -  [(xu  -  l)sin  x2  +  (x2  -  x3)2]  , 
gubject  to  the  constraints: 


0  ^   x1  <  1 

0  <  x2  <  2 

-1  <   x3  ^  1 

1.05  C      x,  <  2 
k 

xx2  +  x22  +  x32  +  xk2     <  3 
x2  +  x  +  x,  +  x2  x.  +  x  x^  +  x^  x^  4.  2 


2  2 

x2   -  x2  Xo  +  x 


1     xh  ~   xl  x32  ^  6' 


For  the  purposes  of  solution  by  Tangent  Search  the  inequalities  given 
above  were  modified  and  expressed  as  eleven  constraint  functions  of  the 
following  form: 

G1(x)  =  xx  >   0 

G2(x)  =  l-x1  >   0 

Go(x)  =  x2  »   0 

Gu(x)  =  2-Xg  2   0 

G   (x)  =  6  -  X2   +  X2X_  -  X2  -  X-^2  Xi  +  X,  Xo    ^   0. 

Klingman  and  Himmelblau  list  the   following  solution  obtained  by  their  method: 
x±  =    .995 
x2  =  U.OU  x  ID"5" 
x3  =   -.953 
xk  =  1.050 
F  =  -U.T93   . 


'Mfe& 


62 


Examination  of  this  result  shows  that  it  lies  very  near  the  intersection 
of  the  3rd,  7th,  9th,  and  10th  constraint  boundaries.   If  it  is  assumed 
that  the  exact  solution  is  at  that  intersection,  the  correct  values  can  be 
immediately  calculated: 


xls  =  +  V .9U875   =   .97^037987 


x2s  =  ° 


[3s  =  -.97^037987... 


x,   =  1.05 

4S 
Fs  =  -4.795  • 

Calculations  with  the  Complex  Method  for  two  initial  random  complexes 
and  with  Tangent  Search  for  three  values  of  r  were  performed  for  Problem 
5,  as  for  all  problems  reported  in  this  section.   When  supplied  with  a 
non-feasible  origin  the  first  phase  of  Tangent  Search  always  correctly  lo- 
cated an  initial  feasible  base  for  the  constrained  minimization  phase. 
This  feature  of  the  algorithm  was  tested  successfully  on  this  and  various 
other  problems;  it  will  not  be  discussed  further. 

All  calculations  reported  in  Tables  Va,  Vb,  and  Vc  began  from  the 
feasible  origin, 

xx  =  .75 

x2  =  -75 

x3  =  -.375 

xh  =   1.3125 

F  =  -1.718179  • 
Initial  step  sizes  for  Tangent  Search  runs  were  .5  for  all  variables. 
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TABLE  Va 
End  Criterion:   F  ^   -U.7 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Const. 

Obj.  Funct. 

Total 

Permissible 

a 

li+719 

10208 

1U669 

10200 

b 

- 

- 

failed 

- 

Tangent  Search 

157 

81 

Base 

Explore 

Alt. 

r  =  .250 

13 

1 

15 

r  =  .375 

279 

125 

15 

1 

21 

r  =  .500 

259 

85 

8 

1 

18 

Table  Vb  provides  terminal  counts  for  the  same  runs.   The  last  two 
lines  give  totals  for  that  point  in  each  solution  at  which  the  base  entered 
a  hypersphere  of  radius  .001  about  the  true  solution.   Both  Complex  Method 
runs  and  the  Tangent  Search  calculation  with  r  =  .25  failed  to  reach 
results  of  comparable  accuracy. 


TABLE  Vb,  Terminal  Results 


Algorithm 

Evaluations 

Trials /Moves 

Note/ 
Case 

Complex 
Meth. 

Const. 

0b  ,1.  Funct. 

Total 

Permissible 

a 

1U879 

10258 

11+827 

10250 

(1) 

b 

10038 

69U2 

9999 

693U 

(2) 

Tangent 
Search 

66h 

303" 

Base 

Explore 

Alt. 

(3) 

r  =  .250 

67 

1 

65 

r  =  .375 

830 

369 

88 

1 

76 

(U) 

r  =  .500 

930 

382 

73 

■k 

8U 

m 

N.  B.  For  final  results  on  failed  cases,  see  Table  Vc  below 

(1)  Failed;  terminated  due  to  excessive  time. 

(2)  Failed;  terminated"  wfyen  trials  reached  9999. 

(3)  Failed:  terminated  when  all  D  less  than  10 
(h)    Succeeded;  results  given  for  the  base  at  wh 


ich  Jt2>;  (x  -  xis)2  -c 

i=l 


001, 


6k 


TABLE  Vc 
Best  Point  Achieved  for  Failed  Runs.   (See  Table  Vb . ) 


Final 
Value 

Case 

(1) 

(2) 

(3) 

xl 

1.000 

.8U0 

.998 

x2 

-.2E-T 

.^03 

1.1E-10 

Xn 

-.9^7 

-.666 

-.9V9 

X), 

1.050 

1.300 

1.050 

F 

-1+.7922 

-2.8903 

-U.7926 

For  this  problem  Tangent  Search  seems  superior  to  both  Multiple- Gradient 
Summation  and  to  the  Complex  Method.   Also,  r  =  .375  seems  the  best  choice 
of  reduction  factor  among  the  three  considered. 
PROBLEM  6 

This  illustration  is  essentially  the  same  as  the  fifth  problem  of 
Klingman  and  Himmelblau.   The  objective  function,  once  again,  originated 
with  Krolak  and  Cooper  (8).   Here  the  minimum  of 

F(w,x,y,z)  =  -  (sin  x  +  z^4^) 
is  desired,  subject  to  constraints, 

0  <  w  <  1 
0  ^-  x  <  2 
0  <-     y     <  1 

0  <   z   <  2 

9     2     2     2 
xd  +  y     +.  z  +  w   <   1  . 

As  before,  for  the  Tangent  Search  solution  the  constraints 

are  re-expressed  as  simple  inequalities  for  vhich  the  boundaries  are 


feasible,  i.e. , 


G  =  w  ZL     0 

G  =  1  -  w  2.  0 
2 


G  =  1  -  x3  -  y2  -  z2  -  w2  2.     0, 
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Klingman  and  Himmelblau  list  the  following  result,  which  is  erroneous: 
w  =  14.89  x  10-6 
x  =  .559 
y  =  U.89  x  10 
z  =  .829 
F  =  -1.530  . 
Examination  of  this  problem  reveals  that  the  constrained  minimum  actual- 
ly lies  at  the  intersection  of  the  1st,  5th,  7th,  and  9th  constraint  bound- 
aries.  That  is,  the  correct  solution  is: 

x  =  1 
s 

w  =  y  =  z   =0 
s   Js    s 

F  =  -1.8U1U71  • •• 
s 

All  calculations  with  both  Tangent  Search  and  Complex  Methods  correctly 
solved  this  problem.   Results  are  given  in  Tables  Via  and  VIb. 

TABLE  Via 
End  Criterion:   F  <  -1.8U 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Constr. 

Ob.i .  Funct . 

Total 

Permissible 

a 

292 

199 

272 

180 

b 

279 

172 

258 

160 

T 

angent  Search 

U31 

222 

Base 

Explore 

Alt. 

r  =  .250 

U3 

3 

37 

r  =  .375 

558 

252 

52 

3 

kQ 

r  =  .500 

56U 

321 

U6 

10 

kQ 
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TABLE  Vlb,  Terminal  Results 


Final 

Algorithm 

Result 
Item 

Complex 

Tangent   Search,    r  = 

a 

b 

.2^0 

.375 

.500 

w 

•98E-9 

1.0E-9 

1.0E-10 

.1+1E-7 

.67E-8 

X 

1.0000 

.99952 

1.0000 

1.0000 

1.0000 

y 

.10E-8 

.98E-9 

1.2E-10 

1.9E-7 

6.9E-8 

z 

•783E-U 

3.11E-2 

1.18E-3 

1.02E-3 

1.15E-3 

F 

-I.8U1U71 

-1.8U1210 

-I.8U1U7I 

-1.81+11+67 

-I.8U1J470 

— — — — — — — 

Evaluations 

————— 
6U0 

7^9 

1U22 

2163 

2016 

Constraints 

Ob.j.Funct. 

1+U0 

1+21 

657 

1060 

1191 

Tr 

ials /Moves 

620 

729 

Total 

Permiss . 

395 

1+07 

- 

- 

- 

Base 

- 

- 

136 

190 

171 

Explore 

- 

- 

5 

22 

39 

Alternate 

- 

- 

119 

160 

1U8 

Note   No. 

(1) 

(2) 

(3) 

(h) 

(h) 

f 


(1)  No  further  significant  changes  in  variables  could  be  found  after  this 
point. 

(2)  Method  could  not  find  a  feasible  replacement  to  the  worst  vertex  at 
this  point. 

(3)  All  step  sizes  become  less  than  10"   at  this  point. 

(1+)  Solution  terminated  due  to  excessive  computing  time  at  this  point; 
final  step  sizes  were  approximately  k   x  10 

The  evidence  of  these  tables  seems  to  show  that  the  Complex  Method  is 
better  than  Tangent  Search  for  solving  this  problem.   The  best  result  for 
the  latter  method  was  with  r  =  .25  .   Finally,  both  of  these  algorithms 
produce  satisfactory  solutions,  even  though  Multiple-Gradient  Summation 
failed. 
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PROBLEM  7 

This   is   the  practical  problem  which  stimulated  development   of  the 
Complex  Method  by  M.    J.    Box    (3).      He  was   unable   to  solve   it  either  by  the 
method  of  Rosenbrock    (T)    or  by  certain   other  techniques. 

The  problem,    as   re-expressed  here,    is   to  minimize  the   following  function 
of   five   independent   variables,   x .  ,    i   =   1,   2,.  ..,5: 


F(x)   =  -  j  j_a2yi  +   a3y2  +  a1+y3  +   a^y^  +   c± 


where, 


+   c-,    -  C2 

-C3(x2  +    .  Olx^)    +  k31  +  k     x2  +  k  oX_ 
+  k3)A  +  k35X5]    xl-2U3U5  +  a^, 


X6  =   ^\  +  k2x2  +  k3x3  +  k^xk  +  k5x5)xx> 

7i  =  k6  +  kTx2  +  k8x3  +  k9xU  +  k10x5    ' 

y2  =  kH  +  k12X2  +  k13X3  +  kl>A  +  kl5X5    »  H 

Y3  =  ki6  +  k17X2  +  kl8x3  +  kl<A  +  k20x5    ' 

^  =  k21  +  k22x2  +  k23X3  +  k2fcxU  +  k25X5    > 

XT  =  (yi  +  y2  +  y3}  Xl  ' 
x8  =  (k26  +  k2Tx2  +  k28x3  +  k29xU  +  k30x5}  Xi  +  x6  +  XT  ; 


subject  to  the  constraints 


0  £  x 
1 

1.2  £  x   £  2.U 

20   £  x3  £  60 

9.0  ±  xh    £  9.3 

6.5  £  x   £  7-0 

o  £  x,  6  29U000 

o 

0  i.       x   £  29U000 
7 

0   -   x_  £  277200   . 
o 
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The  numerical  values  of  a.  ,  ( i=0 ,  1,...,T)   and  k  ,  (j=l,***,  35)   are  to 

1  J 

be    found  in  the   Appendix  of  the  paper  of  M.    J.    Box    (3).      In  the   function 

definition   above,    take 

cn    =   781+0   a, 
1  o 

c     =   100,000   aQ 
c     =   50800   a~    . 
The   initial  point,    specified  by  M.    J.    Box,    is 
xx  =   2.52 
x2  =  2.0 
x3  =  37-5 
xu  =  9.25 
x     =  6.8 

F     =  -2,351,2^3.5, 
which  is    feasible.      According  to  Box  the   correct    solution   is: 

xls   =   k. 537^3 


x2s 

" 

z.n 

X3s 

= 

60 

xHs 

= 

9.3 

X5S< 

= 

7-0 

*Bs 

= 

277, 

,200 

Fs     =  -5,280,33^ 

The  paper  of  Box  (3)  gives  results  achieved  by  him  in  two  trial  runs. 

These  figures  will  be  labelled  "B  ",  and  "B  "  in  the  tables  below.   For 
&  1  '       2 

this  problem  three  initial  random  complexes  were  employed  by  the  present 
author  in  trial  runs  with  the  Complex  Method.   These  will  be  labelled 
a  ,   b  ,  and  c  . 
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A  comparison  of  results  among  these  and  three  Tangent  Search  calculations 
is  given  in  Table  Vila  for  points  comparatively  early  in  each  run.   Step 
sizes  for  the  latter  method  began  at  .1  for  all  variables. 

TABLE  Vila 
Criterion:   Approximately  520  Constraint  Evaluations 


Algorithm 

Evaluations 

Obj. 
Function 

Complex  Method 

Const. 

Obj.  Funct. 

B 

527* 

* 

-5,236,850 

'■  ■  ■■ —  '  1 
a 

538 

U25 

-5,280,169 

b 

527 

U02 

-5,261,088 

c 

521 

U02 

-1^269,^06 

T 

angent  Search 

52U 

261 

-5,273,23^ 

r  =  .250 

r  =  .375 

53k 

2U8 

-5,276,780 

r  =  .500 

518 

2U7 

-5.253, W 

*  At  300  permissible  trials  and  517  total  trials. 

Here  it  seems  that  Tangent  Search  achieves  roughly  comparable  results 
early  in  the  calculation  with  considerably  fewer  objective  function  eval- 
uations.  Comparisons  at  another  convenient  point,  available  for  various 
runs,  are  provided  in  Table  VI lb.  » 

TABLE  Vllb 
Criterion:   Approximately  1200  Constraint  Evaluations 


Algorithm 

Evaluations 

Obj  .Funct. 
(xlO3) 

Note 

Complex  Meth. 

Const. 

Obj.  Funct. 

B~ 

1206 

? 

-5,280.3 

(1) 

c. 

a 

821 

6U5 

-5,280.3 

(2) 

b 

118U 

922 

-5,278.9 

(3) 

c 

1032 

7^5 

-5,271.8 

(3) 

Tan 

gent  Search 

1186 

587 

-5,27^.2 

(U) 

r  =  .250 

r  =  .375 

1201 

538 

-5,280.2 

(5) 

r  =  .500 

1201 

565 

-5,265.1+ 

(6) 

70 


(1)  Reported  "by  Box  as  620  permissible  moves  out  of  1196  trial  moves. 

(2)  Calculation  stopped  after  only  8ll  trials,  U06  of  which  were 
permissible,  with  correct  constrained  minimum! 

(3)  These  calculations  terminated  with  the  results  shown  because  the 
algorithm  was  unable  to  generate  a  feasible  vertex  to  replace  the 
current  worst  vertex.   In  both  cases  the  complex  seemed  to  be  con- 
tracted almost  to  the  limits  of  significance  of  floating-point 
numbers  on  the  CDC  160U.   This  may  indicate  that  larger  over-reflection 
factor  would  be  helpful;  1.3  was  used,  as  recommended. 

(k)    This  is  the  225th  base  point,  found  after  2  conventional  and  83  tangent 
explore  moves. 

(5)  This  is  the  192nd  base  point,  found  after  2  conventional  and  90  tangent 
explore  moves. 

(6)  This  is  the  lU5th  base  point,  found  after  7  conventional  and  90  tangent 
explore  moves. 

Table  VIIc  gives  terminal  results  for  the  Tangent  Search  calculations . 
Included  are  the  results  of  a  run  with  r  =  .375  but  using  a  coupling  cycle 
slightly  different  from  that  described  in  earlier  parts  of  this  paper. 
This  column  is  labelled  "Special".   All  of  these  began  with  initial  step 
sizes  at  .1  for  all  variables.   They  were  terminated  normally. 
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TABLE  VIIc,  Tangent  Search 
End  Criteria:   D.  ^  10  ~°  for  all  i  =  1,  2,...,  n. 


Final 

Normal  Algorithm 

Special 

Result 

r  = 

r  = 
•  375 

.250 

.375 

.500 

xl 

U.6180U 

k. 53810 

U.5UU38 

U.5371^ 

x0 

2.25576 

2.39875 

2.38736 

2.39997 

X 

59.9997 

59.9916 

59.9999 

59.9991 

J 
X), 

9.30000 

9.30000 

9-30000 

9.30000 

Xr 

6.99999 

7.00000 

7.00000 

7.00000 

r 

-527^279 

-5280269 

-52798H 

-5280332 

Moves 

298 

235 

Ul8 

511 

Base 

Explore 

3 

6 

27 

Ik 

Alt. 

107 

12  k 

336 

2kk 

E 

valuations 
Constr. 

1539 

1779 

U88H 

3^23 

Obj.  Funct. 

778 

805 

2092 

1616 

The  special  run  with  r  =  .375  gives  the  best  Tangent  Search  result  for 
this  example.  It  is  almost  exact.  However  the  coupling  cycle  it  uses  was 
much  less  satisfactory  on  other  problems  than  the  one  described  here.  The 
intermediate  step  reduction  size  is  best  with  either  coupling  cycle.  Sat- 
isfactory performance  of  Tangent  Search  seems  marginal  on  this  problem;  it 
is  highly  dependent  on  coupling  modes  and  choice  of  r.  The  Complex  Method 
appears  to  be  somewhat  more  dependable;  however,  it  did  not  work  perfectly 
for  all  initial  random  complexes. 
PROBLEM  8 

This  example  was  constructed  by  the  author  to  test  Tangent  Search's 
special  expedient  for  dealing  with  convex  constraint  boundaries.   The 
objective  surface  is  the  bottom  part  of  the  interior  of  a  hemisphere.   A 
convex  constraint  boundary  lies  between  the  initial  point  and  the  lowest 
point  of  the  surface.   The  constraint  is  chosen  in  such  a  way  that  the  base 
must  always  move  along  it  in  order  to  reach  the  correct  solution. 
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Specifically,  it  is  desired  to  minimize 

F(x,y)  =     y(x-5)2  +  (y-10)2 
subject  to  the  constraint, 

G  (x,y)  =  Ikk  -   (x-5)2  -(y+10)2  Z    0. 
The  initial  point  is 

x  =  11. 5 i  y  =  o,  f  =  11.9269  ... 

and  the  solution  is 

x  =  5,  y  =  2,   F  =  8. 
s        s        s 

Initial  step  sizes  for  Tangent  Search  calculations  were  .1  for  both 
variables.   Results  for  the  usual  Tangent  Search  and  Complex  method  runs 
are  given  in  Table  VIII. 


TABLE  VIII 


End  Criterion 


1:  I/  (x-x  )2  +  (y-y  )2   £  .0001 
V     s         s 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Const. 

Obj.Funct. 

Total 

Permiss . 

a 

138 

61 

131 

60 

b 

1U2 

66 

137 

60 

Tangent  Search 

161 

lU? 

Base 

Expl. 

Alt. 

r  =  .250 

62 

12 

31 

r  =  .375 

1U9 

151 

52 

16 

31 

r  =  .500 

133 

138 

Ul 

17 

29 

The  Complex  Method  seems  more  efficient  for  this  example.   With  Tangent 
Search  the  reduction  factor,  r  =.5,  appears  best. 
PROBLEM  9 

The  Multiple-Gradient  Summation  Method  failed  on  this  example, 
which  corresponds  to  number  seven  of  Klingman  and  Himmelblau  (l).   The 
objective  function  is  the  same  as  that  of  Problem  5  above,  i.e. , 
F  (v,x,y,z)  =  -((z-l)Sin  X  +  (w-y)  ). 
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However,  the  constraints  of  this  problem  are 
0  *  w  <  1 
0  <      x  <   2 

-1  <    y  <  l 

1.05  <  z  <   2 

2    2    2    2 
x  +  y  +  z  +  w      <     2. 

As  before,  these  must  be  re-expressed  in  the  following  manner  for 

Tangent  Search: 

G  =  w  £   0 

1 

G  =  1  -  w  Z      0 
2 

G  =  2-x2  -  y2  -  z2  -  w2  £   0. 
9 

Trial  calculations  with  Tangent   Search   and  the   Complex  Methods  were   initiated 
at  the   feasible   origin, 

w  =  x  =  y  =      .001 

z  =   1.051 

F  =  -.99703  ... 
Initial  step  size  for  the  former  algorithm  was  .0625  for  all  variables. 
Examination  of  the  problem  shows  that  the  constrained  minimum  lies  at  the 
intersection  of  the  3rd,  7th,  and  9th  constraint  boundaries;  specifically, 
this  solution  is: 

ws  =  .66988805... 

xs  =  0 

y  =  -.66988805... 
s 

zs  =  M5 

Fs  =  -2.795. 


lh 


The  Complex  Method  run  beginning  with  random  vertices  chosen  by 
sequence  "a"  failed  to  develop  a  satisfactory  solution.   It  stopped  after 
the  U69th  trial  (262  permissible)  with  the  following  centroid: 
w  =  1.0  x  10-9 
x  =  2.1  x  10"9 
y  =  .9^73648 
z  =  1.050000 
F  =  -1.8975  • • • 
No  significant  changes  in  the  variables  had  taken  place  during  the  last 
eight  permissible  trials;  this  constitutes  the  usual  end  criterion  for  the 
Complex  Method. 

Results  of  the  "b"  calculation  and  the  usual  three  Tangent  Search  runs 
are  given  in  TAble  IXa  for  a  point  closer  to  the  true  solution. 


TABLE  IXa 


End  Criterion: 


(w-w  )   +  (x-x  )2  +  (y-y  )2  +  (z  -  z_)' 
s         s         s  s 


^  .0075 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Constr. 

Obj.Funct. 

Total 

Permiss. 

a 

Failed 

- 

- 

- 

b 

773 

1*28 

751 

1+20 

T 

angent  Search 

503 

263 

Base 

Exp. 

Alt. 

r  =  .250 

100 

1 

uu 

r  =  .375 

U85 

217 

57 

1 

ko 

r  =  .500 

601 

255 

U9 

2 

58 

Ihe  Complex  Method,  run  "b",  did  not  proceed  to  a  more  accurate  solu- 
tion after  the  point  tabulated  above.   Best  results  achieved  for  the  three 
Tangent  Search  runs  are  given  in  Table  IXb. 
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TABLE  IXb,  Tangent  Search 

End  Criteria:   D.  <  10~10  for  all  i  =  1,  2,  ...,  n, 
1 


Fi 

nal 
suit 

Step  Reduction  Factor 

Re 

.250 

.375 

.500 

w 

.67296U 

.670272 

.669901 

X 

5.6E-11 

i+.OE-ll 

1.1E-11 

-.666798 

-.66950U 

-.669875 

z 

1.05000 

1.05000 

1.05000 

F 

-2.79^96 

-2.79500 

-2.79500 

Moves 

169 

170 

103 

Base 

Explore 

1 

2 

6 

Alternate 

98 

135 

lUl 

Evaluations 

10U9 

1^99 

1507 

Constr. 

Ob  j . Funct . 

512 

6^5 

617 

These  results  for  Problem  9  are  clearly  favorable  to  the  Tangent  Search 
Method  over  both  the  Multiple-Gradient  Summation  and  Complex  Methods.   An 
additional  observation  is  that  a  less  accurate  solution  (Table  IXa)  is 
achieved  most  easily  with  r  =  . 375 •   However  the  most  accurate  final 
answer  is  found  (with  fewer  evaluations)  by  r  =  .500. 
PROBLEM  10 

This  illustration  is  taken  from  the  article  of  M.  J.  Box  (3)  where  it 
appears  as  "Problem  B".   As  re-expressed  here,  the  minimum  of 

F  (x)  =  -  (x23  [  9  -  (xx  -  3)2]  /  271TT 
is  required,  under  the  constraints, 
G-l  (x)  =  x±    2   0 


Gg  (x)  =  Xr 


0 


G3  (x)  =  x-j/Ti  -  x2  >   0 
Gh    (x)  =  x3  Z      0 

1 


G^  (x)  =  6  -  X-,  —  0,     where   x„  =  xn  +  T3  x 
j  3    1; 
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The  initial  point,  as  specified  by  Box,  is 

x  =  1,   x  =  .5,     F  =  -.01336  . 
1        2 

Initial  step  sizes  for  Tangent  Search  calculations  were  .1  for  both 
variables.   The  correct  constrained  minimum  is 

x_   =  3,    x   =f37   F  =  -1. 

Is         2s         s 

Each  independent  variable  must  have  explicitly-designated  upper  and 
lower  bounds  for  the  purposes  of  solution  by  the  Complex  Method.   Since 
the  upper  bound  of  x  is  not  stated  by  Box  (3),  the  present  author  chose 
6.0  for  the  value. 

Results  of  the  usual  two  Complex  Method  and  three  Tangent  Search  calcu- 
lations are  given  in  Table  X. 

TABLE  X 
End  Criterion:   F<  -.99995 


Algorithm 

r      Evaluations 

Trials /Moves 

Complex  Meth. 

Constr. 

Obj.Funct. 

Total 

Permiss. 

a 

133 

6k 

122 

60 

b 

122 

.59 

111 

55 

Tangent  Search 

Base 

Expl. 

Alt. 

r  =  .250 

121 

91 

55 

k 

26 

r  =  .375 

153 

110 

51 

7 

36 

r  =  .500 

159 

107 

38 

10 

36 

The  Complex  Method  seems  to  solve  this  problem  more  economically. 
With  Tangent  Search  the  smallest  step  reduction  factor  seems  the  best 
choice. 
PROBLEM  11 

This  is  the  Post  Office  Parcel  Problem,  originally  proposed  by 
Rosenbrock  (7)>  but  here  subject  to  the  constraints  imposed  by  Box(3)< 

Expressed  as  a  minimization  problem,  the  objective  function  is 
F(x)  =  -x-j^  x2  x   , 
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and  the   constraints   to  be   observed  are, 
0     6r     v       £     20 
0    £     x2    6.     11 

0      ^     x        ^     142 

~       3     - 

xn    +  2(x     +  x^»)   £    72    . 
1  2  ° 

The  designated  origin  for  all  computations  is 

x,  =  18,   x  =10,   x  =16,   F  =  -2880  . 
1      '    2         3 

Initial  step  size  on  Tangent  Searches  was  unity  for  all  variables. 

The  correct  answer  for  Problem  11  is 

x   =20,   x   =  11,   x   =15,   F  =  -3300. 
Is  2s  3S  s 

The  usual  trial  runs  with  Complex  and  Tangent  Search  algorithms  were 
made.   A  comparison  of  results  is  provided  by  Table  XI.   As  stated  prev- 
iously ?  the  centroid  (or  corresponding  objective  function  value)  is  the 
reference  compared  with  the  end  criterion  for  Complex  Method  calculations 
The  base  point  is  similarly  the  reference  for  Tangent  Search. 

TABLE  XI 
End  Criterion:   F  *.  -3299-9 


Algorithm 

Evaluations 

Trials /Moves 

Complex  Meth. 

Constr. 

Obj .Funct. 

Total 

Permiss . 

a 

289 

209 

271 

150 

b 

368 

258 

359 

190 

Tangent  Search 

195 

107 

Base 

Expl. 

Alt. 

r  =  .250 

38 

1 

30 

. r  =  .375 

186 

97 

25 

1 

28 

r  =  .500 

268 

131 

27 

3 

U0 

The  evidence  of  Table  XI  seems  favorable  to  Tangent  Search,  and  with 
that  method  the  intermediate  step  reduction  factor  produces  the  best  result, 
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Conclusions 

The  following  general  summary  of  the  results  of  trial  computations 
tabulated  above  is  based  to  some  degree  upon  qualitative  judgments.    Results 
with  different  methods  for  a  particular  problem  are  called  "similar"  if  the 
total  number  of  function  evaluations  required  to  reach  a  valid  solution  of 
pre-defined  accuracy  are  approximately  the  same.   This  total  for  each  run 
includes  both  objective  and  constraint  function  calls. 

Otherwise,  a  method  is  termed  "superior"  to  another  for  anv  particular 
problem  if  a  valid  solution  is  attained  with  substantially  fewer  evalu- 
ations.  Also,  of  course,  a  method  is  superior  if  it  achieves  a  valid  con- 
strained minimum  while  the  compared  method  fails  to  do  so.   Finally, 
Tangent  Search  is  never  rated  superior  to  another  if  a  better  result  is 
achieved  with  only  one  of  the  step  reduction  factors  considered.   The  con- 
verse concepts  define  when  one  algorithm  is  considered  inferior  to  another. 
Table  XII  provides  a  summary  of  these  judgments  applied  to  the  tables  of 

the  previous  section. 

TABLE  XII 

Number  of  Examples  for  Which  Tangent  Search  Was  Similar,  Superior,  or 

Inferior  to  Other  Methods 


Qualitative 
Comparison 

Multiple  Gradient 
Summat  i  on 

Sequential 
Search 

Complex 
Method 

Similar 

1 

1 

1 

Superior 

k 

1 

6 

Inferior 

0 

0 

h 

Total  Cases 

5 

2 

11 
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It  appears  from  this  summary  that  the  Tangent  Search  Method 
can  be  considered  a  reasonably-efficient  means  of  locating  a  constrained 
minimum,  provided  the  other  methods  are  so  classified.  However,  it  must 
be  remarked  that  all  but  one  of  the  examples  considered  are  text-book  or 
artificial  problems.  Problem  7  is  the  only  practical  example  and  for  it 
the  performance  of  the  Complex  Method  was  somewhat  superior. 

Another  general  conclusion  would  be  that  with  Tangent  Search  a  step 
reduction  factor  of  .25  or  .375  seems  better  as  a  generally-applicable 
choice  than  . 5  •   On  only  one  problem,  the  8th,  was  the  largest  factor 
superior  in  performance  to  both  of  the  others.   It  also  appears  that 
calculations  with  r  =  .25  converged  toward  the  solution  more  rapidly  than 
the  others  at  the  beginning  of  problems  in  k   or  5  dimensions.   However, 
more  accurate  final  results  were  often  obtained  for  runs  with  r  =  . 375  or 
.500  (at  the  cost  of  extra  evaluations).   In  particular,  the  best  final 
values  for  Problems  5  and  7  were  achieved  with  r  =  .  3755  and  best  final 
values  for  Problems  6  and  9  were  for  r  =  .500  .   These  examples  of  higher 
dimensionality  were  the  most  troublesome  considered.   The  step  reduction 
factor  r  =  .25  seemed  to  work  well  for  problems  of  dimensionality  3  or  less. 

As  a  final  conclusion,  based  upon  intuition  alone,  the  author  wants  to 
endorse  an  often-stated  opinion  of  most  other  investigators  in  this  area; 
viz. ,  the  final  word  on  the  best  method  of  obtaining  constrained  minima/ 
optima  remains  to  be  said.   The  method  of  Davidon  (9)»  the  foundation  of 
the  method  of  Fletcher  and  Powell  (10),  has  emerged  as  a  tentative  "optimum 
optimizer"  for  the  unconstrained  problem.   However,  no  such  leading  can- 
didate for  solving  the  general  problem  of  non-linear  mathematical  program- 
ming has  yet  been  heralded. 
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Future  Work 

Several  avenues  for  further  investigations  have  been  opened  in  the 
previous  discussion.   The  present  author  hopes  to  investigate  some  of  these: 

(1)  The  step  reduction  factor  in  Tangent  Search  probably  should  be  sub- 
ject to  modification  by  the  algorithm,  itself.   Probably  it  should 
start  at  the  smaller  end  of  some  range  and  become  larger  in  the  vicinity 
of  a  solution. 

(2)  Preliminary  results  of  computational  experiments  now  in  progress  show 
that  partial  derivatives  should  be  calculated  by  the  method  of  Wengert(5) 
or,  perhaps,  that  of  Smith  (6).   A  package  of  subroutines  allowing  users 
to  utilize  the  former  method  conveniently  has  been  developed  at  the 
Computer  Facility,  U.  S.  Naval  Postgraduate  School.   All  necessary  jumps 
and  calling  lists  are  generated  by  the  Fortran-63  Compiler;  therefore, 
the  user  is  not  inconvenienced  by  the  fact  that  the  derivatives  are  ac- 
tually developed  in  subroutines.   Of  course  the  constraint  functions  must 
be  capable  of  explicit  statement  if  such  a  method  is  to  be  used. 

(3)  Explicit   constraints  should  be  treated  separately  to  avoid  useless 
evaluation  of  partial  derivatives  which  are  always  either  zero  or  unity. 
For  the  same  reason  the  procedure  which  calculates  derivatives  of  the 
implicit  constraints  could  utilize  a  simple  switch  so  that  only  those 
presently  needed  would  be  evaluated. 

(h)      More  computational  comparisons  between  Sequential  Search  and  Tangent 
Search  will  be  necessary  to  establish  whether  or  not  one  method  is  gen- 
erally superior. 

(5)   The  author  feels  that  more  investigation  of  the  parameters  oc.  and  k  is 
necessary  in  establishing  their  best  general  values  in  the  Complex  Method. 
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M.  J.  Box  chose  °C  =  1.3  and  k  =  2n  on  the  basis  of  results  after  only  200 
trials  on  only  two  examples.   However,  in  several  of  the  calculations  re- 
ported above  a  complex  became  totally  contracted  before  the  correct  solution 
was  attained.   This  might  be  remedied  by  using  a  larger  value  for  the  over- 
reflection  factor  or  more  vertices. 
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